Application Project - New York City Bike Sharing
New York City Bike Project
Überblick
In diesem Anwendungsprojekt werden Daten der City of New York bearbeitet, visualisiert und analysiert, die Informationen zum städtischen Bikesharing aus dem Jahr 2016 enthalten. Das Hauptaugenmerk in diesem Projekt liegt darauf möglichst sinnvoll und genau die Anzahl der Nutzer an den verschiedenen Stationen vorherzusagen. Zu diesem Zweck werden zusätzlich Wetterdaten der Stadt New York hinzugezogen.
Die Vorhersage der Nutzerzahlen geschieht mit Hilfe eines künstlichen neuronalen Netzes, einem sogenannten Long Short Term Memory Netz, kurz LSTM. Um die Performance des Netzes zu überprüfen wird das Ergebnis mit dem Ergebnis einer linearen Regression nach der Methode der kleinsten Quadrate, auch OLS Regression genannt, verglichen.
Daten
Die Daten zum Bikesharing sind folgendermaßen strukturiert. Jede Zeile des ursprünglichen Datensatzes repräsentiert eine Fahrt mit einem Fahrrad der Bikesharing Flotte. Zu jeder Fahrt ist die Dauer der Leihe, der Start- und Endzeitpunkt, sowie die Start- und Endstation festgehalten. Außerdem ist die Indentifikationsnummer des jeweiligen Fahrrades vorhanden. Des Weiteren sind Kundeninformationen festgehalten. Das Geschlecht und Geburtsjahr, sowie der Nutzertyp. Bei Letzterem wird zwischen Customer und Subscriber unterschieden. Subscriber haben ein jährliches Abonnement abgeschlossen, wohingegen Customer einen Ein- oder Drei-Tages-Pass erworben haben. Zu Analysezwecken wird dieser ursprüngliche Datensatz im Laufe des Projektes aggregiert. Dies geschieht auf verschiedenen Ebenen. Zum Einen werden die Daten nach Stunde und Tag und zum Anderen nach Leihstation und Stadtbezirk zusammengefasst.
Die täglichen Wetterdaten halten die täglichen Maximal- und Minimaltemperaturen fest, sowie die Tagesdurchschnittstemperatur. Des Weiteren geben Sie Auskunft über den gefallenen Regen, sowie den neuen Schneefall und die Tiefe der aktuellen Schneedecke.
Die stündlichen Wetterdaten geben Auskunft über die Temperatur, die Windgeschwindigkeit, sowie den gefallenen Regen.
Forschungsfragen
Nach einigen Diskussionen und Analysen der Daten haben sich zwei interessante Forschungsfragen ergeben.
Können die Nutzerzahlen pro Station und Stunde zuverlässig vorausgesagt werden? und
Können die Nutzerzahlen pro Station und Tag zuverlässig vorausgesagt werden?
Im späteren Verlauf des Projektes haben sich aufgrund der Datenmenge Probleme in den Analysen ergeben. Aufgrunddessen werden die Forschungsfragen wie folgt bearbeitet:
1. Können die Nutzerzahlen pro Bezirk und Stunde zuverlässig vorausgesagt werden? und
2. Können die Nutzerzahlen pro Bezirk und Tag zuverlässig vorausgesagt werden?
Durch die Aggregation auf Bezirksebene verringert sich die Anzahl an Variablen erheblich. Aufgrund der dennoch enorm großen Datenmengen und den daraus resultierenden Problemen, werden die Vorhersagen mit dem LSTM Netz und der linearen Regression exemplarisch für einen Bezirk durchgeführt.
Preprocessing
In einem ersten Schritt, werden alle zunächst zur Verfügung stehenden Datensätze geladen und entsprechend der Analysezwecke aufbereitet.
Zunächst werden die Datensätze geladen, die die Wetterdaten enthalten. Je Forschungsfrage gibt es einen eigenen Datensatz. Der Datensatz weather enthält tägliche Wetterdaten, wohingegen der Datensatz hourly_weather stündliche Wetterdaten enthält.
weather <- read.csv("Data/weather_data_nyc_centralpark_2016(1).csv")
hourly_weather <- read.csv("Data/hourly_weather.csv")Da die Bikesharingdaten von der Stadt New York monatlich gespeichert werden, werden sie zunächst je Monat in einen Data Frame geladen.
januar <- read_csv("H:/Projekt Daten/201608-citibike-tripdata/201601-citibike-tripdata.csv")
februar <- read_csv("H:/Projekt Daten/201608-citibike-tripdata/201602-citibike-tripdata.csv")
maerz <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201603-citibike-tripdata.csv")
april <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201604-citibike-tripdata.csv")
mai <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201605-citibike-tripdata.csv")
juni <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201606-citibike-tripdata.csv")
juli <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201607-citibike-tripdata.csv")
august <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201608-citibike-tripdata.csv")
september <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201609-citibike-tripdata.csv")
oktober <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201610-citibike-tripdata.csv")
november <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201611-citibike-tripdata.csv")
dezember <-read_csv("H:/Projekt Daten/201608-citibike-tripdata/201612-citibike-tripdata.csv")Da es in R schwierig sein kann mit Spaltennamen zu arbeiten, die Leerzeichen enthalten, werden diese durch Unterstriche ersetzt.
names(januar) <- str_replace_all(names(januar), c(" " = "_"))
names(februar) <- str_replace_all(names(februar), c(" " = "_"))
names(maerz) <- str_replace_all(names(maerz), c(" " = "_"))
names(april) <- str_replace_all(names(april), c(" " = "_"))
names(mai) <- str_replace_all(names(mai), c(" " = "_"))
names(juni) <- str_replace_all(names(juni), c(" " = "_"))
names(juli) <- str_replace_all(names(juli), c(" " = "_"))
names(august) <- str_replace_all(names(august), c(" " = "_"))
names(september) <- str_replace_all(names(september), c(" " = "_"))
names(oktober) <- str_replace_all(names(oktober), c(" " = "_"))
names(november) <- str_replace_all(names(november), c(" " = "_"))
names(dezember) <- str_replace_all(names(dezember), c(" " = "_"))In einem nächsten Schritt werden die Spaltennamen der Monate Oktober bis Dezember an die Namen der anderen Monate angepasst, damit ein Datensatz entstehen kann, der alle Monate umfasst.
dezember <- dezember %>%
rename(tripduration = Trip_Duration,
starttime = Start_Time,
stoptime = Stop_Time,
start_station_id = Start_Station_ID,
start_station_name = Start_Station_Name,
start_station_latitude = Start_Station_Latitude,
start_station_longitude = Start_Station_Longitude,
end_station_id = End_Station_ID,
end_station_name = End_Station_Name,
end_station_latitude = End_Station_Latitude,
end_station_longitude = End_Station_Longitude,
bikeid = Bike_ID,
usertype = User_Type,
birth_year = Birth_Year,
gender = Gender)
november <- november %>%
rename(tripduration = Trip_Duration,
starttime = Start_Time,
stoptime = Stop_Time,
start_station_id = Start_Station_ID,
start_station_name = Start_Station_Name,
start_station_latitude = Start_Station_Latitude,
start_station_longitude = Start_Station_Longitude,
end_station_id = End_Station_ID,
end_station_name = End_Station_Name,
end_station_latitude = End_Station_Latitude,
end_station_longitude = End_Station_Longitude,
bikeid = Bike_ID,
usertype = User_Type,
birth_year = Birth_Year,
gender = Gender)
oktober <- oktober %>%
rename(tripduration = Trip_Duration,
starttime = Start_Time,
stoptime = Stop_Time,
start_station_id = Start_Station_ID,
start_station_name = Start_Station_Name,
start_station_latitude = Start_Station_Latitude,
start_station_longitude = Start_Station_Longitude,
end_station_id = End_Station_ID,
end_station_name = End_Station_Name,
end_station_latitude = End_Station_Latitude,
end_station_longitude = End_Station_Longitude,
bikeid = Bike_ID,
usertype = User_Type,
birth_year = Birth_Year,
gender = Gender)Aufgrund der Datenmenge werden die Monate einzeln abgespeichert und nach dem Laden zu einem Data Frame zusammengefasst. So ist während der Projektarbeit sichergestellt, dass alle Teammitglieder zu jeder Zeit über GitHub die jeweils aktuellen Daten nutzen können und keine Diskrepanzen entstehen.
saveRDS(januar, "Bike_Data_Jan.RDS")
saveRDS(februar, "Bike_Data_Feb.RDS")
saveRDS(maerz, "Bike_Data_Mrz.RDS")
saveRDS(april, "Bike_Data_Apr.RDS")
saveRDS(mai, "Bike_Data_Mai.RDS")
saveRDS(juni, "Bike_Data_Jun.RDS")
saveRDS(juli, "Bike_Data_Jul.RDS")
saveRDS(august, "Bike_Data_Aug.RDS")
saveRDS(september, "Bike_Data_Sep.RDS")
saveRDS(oktober, "Bike_Data_Okt.RDS")
saveRDS(november, "Bike_Data_Nov.RDS")
saveRDS(dezember, "Bike_Data_Dez.RDS")Die einzelnen Datensätze für jeden Monat werden geladen und in einem Datensatz mit dem Namen “bike” zusammengefasst. Da in den Datensätzen Januar bis September die Start- und Stopzeiten als character formatiert sind, in den Monaten Oktober bis Dezember aber als Unixtime bzw. bereits richtig umgewandelt als datetime, muss vor der Zusammenführung das Format vereinheitlicht werden. Hierfür wird in den Monaten Januar bis September der Datentyp von character zu datetime umgewandelt. So ist weiteres Arbeiten mit den Datums- und Zeitangaben problemlos möglich.
df1 <- read_rds("Data/Bike_Data_Jan.RDS")
df1$starttime <- mdy_hms(df1$starttime)
df1$stoptime <- mdy_hms(df1$stoptime)
df2 <- read_rds("Data/Bike_Data_Feb.RDS")
df2$starttime <- mdy_hms(df2$starttime)
df2$stoptime <- mdy_hms(df2$stoptime)
df3 <- read_rds("Data/Bike_Data_Mrz.RDS")
df3$starttime <- mdy_hms(df3$starttime)
df3$stoptime <- mdy_hms(df3$stoptime)
df4 <- read_rds("Data/Bike_Data_Apr.RDS")
df4$starttime <- mdy_hms(df4$starttime)
df4$stoptime <- mdy_hms(df4$stoptime)
df5 <- read_rds("Data/Bike_Data_Mai.RDS")
df5$starttime <- mdy_hms(df5$starttime)
df5$stoptime <- mdy_hms(df5$stoptime)
df6 <- read_rds("Data/Bike_Data_Jun.RDS")
df6$starttime <- mdy_hms(df6$starttime)
df6$stoptime <- mdy_hms(df6$stoptime)
df7 <- read_rds("Data/Bike_Data_Jul.RDS")
df7$starttime <- mdy_hms(df7$starttime)
df7$stoptime <- mdy_hms(df7$stoptime)
df8 <- read_rds("Data/Bike_Data_Aug.RDS")
df8$starttime <- mdy_hms(df8$starttime)
df8$stoptime <- mdy_hms(df8$stoptime)
df9 <- read_rds("Data/Bike_Data_Sep.RDS")
df9$starttime <- mdy_hms(df9$starttime)
df9$stoptime <- mdy_hms(df9$stoptime)
df10 <- read_rds("Data/Bike_Data_Okt.RDS")
df11 <- read_rds("Data/Bike_Data_Nov.RDS")
df12 <- read_rds("Data/Bike_Data_Dez.RDS")
bike <- rbind(df1, df2, df3, df4, df5, df6, df7, df8, df9, df10, df11, df12)
rm(df1, df2, df3, df4, df5, df6, df7, df8, df9, df10, df11, df12)Missing Values
Im nächsten Schritt des Preprocessing wird in den Datensätzen bike und weather nach fehlenden Werten gesucht und diese entfernt oder ersetzt. Dies ist notwendig, da die späteren Analyseergebnisse durch fehlende Werte beeinträchtigt werden können. Die fehlenden Werte werden pro Spalte aufsummiert.
Im Datensatz der täglichen Wetterdaten befindet sich in den Spalten Niederschlagsmenge, Schneefall und Schneetiefe anstelle von numerischen Werten and einigen Stellen der Wert T. Dieser steht für eine nicht messbare Menge an Regen oder Schnee. Würde der Wert ein Platzhalter für fehlende Werte darstellen, könnten diese wie unten vorgestellt mit Hilfe der Spline-Interpolation sinnvoll geschätzt werden. Diese Variante wird nicht weiter genutzt, da der Wert T keinen fehlenden Wert darstellt.
# Überprüfe auf missing values
apply(bike, 2, function(x) sum(is.na(x)))
apply(weather, 2, function(x) sum(is.na(x)))
# Mögliche Variante für missing values beim Wetter
weather_test <- weather
is.na(weather_test$precipitation) <- which(weather$precipitation == 'T')
is.na(weather_test$snow.fall) <- which(weather$snow.fall == 'T')
is.na(weather_test$snow.depth) <- which(weather$snow.depth == 'T')
weather_test$precipitation <- as.numeric(weather_test$precipitation)
weather_test[is.na(weather_test$precipitation),"precipitation"] <- 0.01
weather_test_interpolated <- na_interpolation(weather_test$precipitation, option = "spline")Die Spalten usertype und birth_year des bike Datensatzes enthalten als einzige fehlende Werte. Um einen sauberen Datensatz zu erhalten, werden die Zeilen, die fehlende Werte enthalten aus dem Datensatz entfernt.
Die hat zwei Gründe. Zum Einen kann das Geburtsjahr, sowie auch der Nutzertyp (Subscriber oder Customer) nicht sinnvoll interpoliert werden. Es können von den davor und danach liegenden Datenpunkten keine sinnvollen Rückschlüsse auf die fehlenden Werte gezogen werden. Zum Anderen liegt hier ein ausreichend großer Datensatz vor (über 13 000 000 Zeilen), sodass der Informationsverlust, der durch das Löschen von ca. 1 600 000 Zeilen gering ist. So können Verzerrungen durch falsch interpolierte Werte vermieden werden.
Aufbereitung
Tägliche Wetterdaten
Um mit den täglichen Wetterdaten arbeiten zu können, muss zunächst die Spalte date in das richtige Datumsformat überführt werden. Die Temperaturangaben in Fahrenheit werden in Celsius, die Menge des Regens und des Schnees in Inches in Millimeter umgerechnet. Kann die Menge des gefallenen Regens oder Schnees nicht gemessen werden, wird dies durch ein T gekennzeichnet. Um mit den Daten problemlos arbeiten zu können, wird das T durch den Wert 0.01 ersetzt. So ist ein numerischer Wert vorhanden, der kennzeichnet, dass Regen oder Schnee gefallen ist. Jedoch ist dieser sehr gering, sodass er keinen unverhältnismäßigen Einfluss auf etwaige Analyseergebnisse hat. Die T-Werte werden zunächst durch den Wert 100 ersetzt, da so Probleme beim Umrechnen umgangen werden, die durch z.B. NAs entstehen. In keiner der Spalten kommt vorher der Wert 100 vor, sodass nach dem Umrechnen weiterhin klar ist, welche Werte vorher ein T waren. Nach dem Umrechnen haben diese den Wert 2540 und können problemlos mit 0.01 ersetzt werden.
# Datum richtig formatieren
weather$date <- dmy(weather$date)
# Temperatur in Celsius umrechnen
weather <- weather %>%
mutate(maximum.temperature = round((maximum.temperature-32)*5/9, 2),
minimum.temperature = round((minimum.temperature-32)*5/9, 2),
average.temperature = round((average.temperature-32)*5/9, 2))
levels(weather$precipitation)
levels(weather$snow.fall)
levels(weather$snow.depth)
# Niederschlagsmenge formatieren, T-Werte durch 100 ersetzen und zurückformatieren
weather$precipitation <- as.character(weather$precipitation)
weather$precipitation[weather$precipitation == "T"] <- 100
weather$precipitation <- as.double(weather$precipitation)
# Schneefall und -tiefe formatieren, T-Werte durch 100 ersetzten und zurückformatieren
weather$snow.fall <- as.character(weather$snow.fall)
weather$snow.fall[weather$snow.fall == "T"] <- 100
weather$snow.fall <- as.double(weather$snow.fall)
weather$snow.depth <- as.character(weather$snow.depth)
weather$snow.depth[weather$snow.depth == "T"] <- 100
weather$snow.depth <- as.double(weather$snow.depth)
# Regen und Schnee in mm umrechen
weather <- weather %>%
mutate(precipitation = round(precipitation*25.4, 2),
snow.fall = round(snow.fall*25.4, 2),
snow.depth = round(snow.depth*25.4, 2))
# Die ehemaligen T-Werte (jetzt 2540) mit 0.01 ersetzten
max(weather$precipitation)
weather$precipitation[weather$precipitation == 2540] <- 0.01
max(weather$snow.fall)
weather$snow.fall[weather$snow.fall == 2540] <- 0.01
max(weather$snow.depth)
weather$snow.depth[weather$snow.depth == 2540] <- 0.01
# Wetterdaten speichern
saveRDS(weather, "Data/weather_daily_2016.rds")Die Wetterdaten für das Jahr 2017 werden genauso aufbereitet wie die Vorherigen. Die Maßeinheiten Fahrenheit und Inches werden wieder in Celsius und Millimeter umgerechnet. Da die Spalte mit der täglichen Durchschnittstemperatur ausschließlich fehlende Werte enthält, wird sie aus dem Mittelwert der minimal und maximal Temperatur ermittelt. Dies ist ebenfalls das Vorgehen des National Climatic Data Center, kurz NCDC, woher die Daten stammen. Die Spalten Niederschlagsmenge, Schneefall und Schneetiefe enthalten im Jahr 2017 keine T-Werte wie noch im Jahr 2016 und müssen hier nicht dementsprechend aufbereitet werden.
# Daten laden
weather17 <- read_csv("Data/daily_weather_2017.csv")
# Daten für 2017 auswählen, Spaltennamen dem Vorjahr anpassen und Durchschnittstemperatur berechnen
weather17 <- weather17 %>%
filter(year(DATE) == 2017) %>%
select(date = DATE,
maximum.temperature = TMAX,
minimum.temperature = TMIN,
precipitation = PRCP,
snow.fall = SNOW,
snow.depth = SNWD) %>%
mutate(average.temperature = (minimum.temperature+maximum.temperature)/2)
# Fehlende Werte überprüfen
apply(df, 2, function(x) sum(is.na(x)))
# Einheiten umrechnen
weather17 <- weather17 %>%
mutate(maximum.temperature = round((maximum.temperature-32)*5/9, 2),
minimum.temperature = round((minimum.temperature-32)*5/9, 2),
average.temperature = round((average.temperature-32)*5/9, 2),
precipitation = round(precipitation*25.4, 2),
snow.fall = round(snow.fall*25.4, 2),
snow.depth = round(snow.depth*25.4, 2))
# Data Frame wie Vorjahr sortieren
weather17 <- weather17[,c(1:3,7,4:6)]
# Datensatz speichern
saveRDS(weather17, "Data/weather_daily_2017.RDS")Stündliche Wetterdaten
Auch die stündlichen Wetterdaten müssen vor Gebrauch bearbeitet werden. Der Datensatz hat in vielen Spalten sehr viele fehlende Werte. In einem ersten Schritt werden die Spalten ausgesucht, die weiterhin betrachtet werden sollen. Die Spalten Datum, Temperatur in Grad Celsius, Windgeschwindigkeit in Km/h, sowie die binären Spalten Nebel, Regen, Schnee, Hagel, Gewitter und Tornado werden weiterhin genutzt. Die Spalte, die z.B den gefallenen Regen in mm enthält, weißt mehr fehlende Werte auf, als Werte da sind. Da dadurch ein sinnvolles interpolieren der Werte nicht mehr möglich ist, wird diese Variable, sowie einige andere von der weiteren Betrachtung ausgeschlossen. Die Spalten Temperatur und Windgeschwindigkeit weisen einige fehlende Werte auf, die aber aufgrund ihrer geringen Anzahl sinnvoll interpoiert werden können und werden. Hierfür wird die Spline-Interpolation genutzt. In diesem Datensatz wird die Spalte, die Datum und Uhrzeit enthält aufgespalten, sodass Datum und Stunde als getrennte Spalten übrig bleiben. Im letzten Schritt der Aufbereitung werden doppelte Beobachtungen behandelt. Für manche Stunden sind verschiedene Temperaturen und Windgeschwindigkeiten bekannt, die sich jeweils nur marginal unterscheiden. Es wird jeweils der erste gemessene Wert für die weiteren Analysen beibehalten.
# Nur Wetterdaten aus 2016 behalten
hourly_weather <- hourly_weather %>%
filter(date(pickup_datetime) > "2015-12-31" & date(pickup_datetime) < "2017-01-01")
# Fehlende Werte überprüfen
apply(hourly_weather, 2, function(x) sum(is.na(x)))
# Spalten auswählen
hourly_weather <- hourly_weather %>%
select(pickup_datetime, tempm, wspdm, fog, rain, snow, hail, thunder, tornado)
# Auf fehlende Werte überprüfen und diese interpolieren
apply(hourly_weather, 2, function(x) sum(is.na(x)))
hourly_weather[is.na(hourly_weather$tempm),]
hourly_weather$tempm <- na_interpolation(hourly_weather$tempm, option = "spline")
hourly_weather[is.na(hourly_weather$wspdm),]
hourly_weather$wspdm <- na_interpolation(hourly_weather$wspdm, option = "spline")
apply(hourly_weather, 2, function(x) sum(is.na(x)))
# Spalten umbenennen
hourly_weather <- hourly_weather %>%
rename(datetime = pickup_datetime,
temp = tempm,
wspd = wspdm)
# Datum und Stunde trennen
hourly_weather <- hourly_weather %>%
mutate(date = as_date(datetime),
hour = hour(datetime))
hourly_weather$datetime <- NULL
# Reihenfolge der Spalten festlegen und nur einen Messwert je Stunde behalten
hourly_weather <- hourly_weather[,c(9,10,1:8)]
hourly_weather <- hourly_weather %>% distinct(date, hour, .keep_all = T)
# Daten speichern
saveRDS(hourly_weather, "hourly_weather.RDS")Aufgrund der dünnen Datenlage zum stündlichen Wetter im Jahr 2017, werden zwei andere Datensätze für die Jahre 2016 und 2017 genutzt, die das stündliche Wetter enthalten.
Auch diese müssen zunächst aufbereitet werden. Da sehr viele Variablen vorhanden sind, die zum Teil viele fehlende Werte enthalten werden die stündlichen Daten für die Temperatur an der Messstation in Fahrenheit, die Windgeschwindigkeit in Meilen pro Stunde, sowie die Niederschlagsmenge in Inches ausgewählt. Fahrenheit wird in Celsius, Meilen pro Stunde in Kilometer pro Stunde und Inches in Millimeter umgerechnet. Kann die Menge des gefallenen Regens nicht gemessen werden, wird dies durch ein T gekennzeichnet. Um mit den Daten problemlos arbeiten zu können, wird das T durch den Wert 0.01 ersetzt. So ist ein numerischer Wert vorhanden, der kennzeichnet, dass Regen gefallen ist. Jedoch ist dieser sehr gering, sodass er keinen unverhältnismäßigen Einfluss auf etwaige Analyseergebnisse hat. Die T-Werte werden zunächst durch den Wert 1000 ersetzt, da so Probleme beim Umrechnen umgangen werden, die durch z.B. NAs entstehen. In keiner der Spalten kommt vorher der Wert 1000 vor, sodass nach dem Umrechnen weiterhin klar ist, welche Werte vorher ein T waren. Nach dem Umrechnen haben diese den Wert 25400 und können problemlos mit 0.01 ersetzt werden. Sind mehrere Werte für eine Stunde vorhanden, so wird nur der erste Wert weiter berücksichtigt. Verbliebene fehlende Werte werden mit Hilfe der Spline-Interpolation ersetzt. Jeder Tag endet mit einem Report um 23:59:00 Uhr, der ausschließlich fehlende Werte enthält. Diese Zeilen werden aus dem Datensatz entfernt. Dies geschieht für die Jahre 2016 und 2017 getrennt.
# Daten einlesen
df <- read_csv("Data/hourly_weather_new.csv")
# stündliche Daten auswählen
df <- df %>% select(DATE, HourlyAltimeterSetting, HourlyDewPointTemperature, HourlyDryBulbTemperature,
HourlyPrecipitation, HourlyPresentWeatherType, HourlyPressureChange, HourlyPressureTendency,
HourlyRelativeHumidity, HourlySeaLevelPressure, HourlySkyConditions, HourlyStationPressure,
HourlyVisibility, HourlyWetBulbTemperature, HourlyWindDirection, HourlyWindGustSpeed,
HourlyWindSpeed)
# Fehlende Werte überprüfen
apply(df, 2, function(x) sum(is.na(x)))
# Temperatur, Windgeschwindigkeit und Niederschlagsmenge auswählen
df <- df %>% select(DATE, HourlyDryBulbTemperature, HourlyPrecipitation, HourlyWindSpeed)
# Datensatz nach Jahren 2016 und 2017 aufsplitten
df16 <- df %>%
filter(year(DATE)==2016)
df17 <- df %>%
filter(year(DATE)==2017)
# Jeweils fehlende Werte überprüfen
apply(df16, 2, function(x) sum(is.na(x)))
apply(df17, 2, function(x) sum(is.na(x)))
# 2016
# Datum und Stunde trennen, sowie Uhrzeit als eigene Spalte behalten
df16 <- df16 %>%
mutate(date = as_date(DATE),
hour = hour(DATE),
time = paste(hour(DATE), minute(DATE), second(DATE), sep = ":"))
# Alle Zeilen mit Uhrzeit 23:59:00 entfernen, da letzter Tagesreport mit allen Werten NA
df16 <- df16 %>%
filter(time != "23:59:0")
df16$time <- NULL
df16$DATE <- NULL
# Temperatur und Windgeschwindigkeit umrechnen
df16 <- df16 %>%
mutate(temp = round((HourlyDryBulbTemperature-32)*5/9, 2),
wdsp = round(HourlyWindSpeed*1.609344, 2))
df16$HourlyDryBulbTemperature <- NULL
df16$HourlyWindSpeed <- NULL
#Niederschlagsmenge umbenennen, Werte mit Buchstaben korrigieren, T durch 1000 ersetzten, in mm umrechnen und die ehemaligen T-Werte in 0.01mm ändern
df16 <- df16 %>%
rename(precip = HourlyPrecipitation)
df16$precip <- str_replace_all(df16$precip, c("s" = ""))
df16$precip[df16$precip == "T"] <- 1000
df16$precip <- as.double(df16$precip)
df16 <- df16 %>%
mutate(precip = round(precip*25.4, 2))
max(df16$precip, na.rm = T)
df16$precip[df16$precip == 25400] <- 0.01
# Reihenfolge der Spalten ändern und nur eine Zeile pro Stunde behalten
df16 <- df16[,c(2,3,4,5,1)]
df16 <- df16 %>%
distinct(date, hour, .keep_all = T)
# Fehlende Werte interpolieren
df16[,3:5] <- na_interpolation(df16[,3:5], option = "spline")
apply(df16, 2, function(x) sum(is.na(x)))
# 2017
# Datum und Stunde trennen, sowie Uhrzeit als eigene Spalte behalten
df17 <- df17 %>%
mutate(date = as_date(DATE),
hour = hour(DATE),
time = paste(hour(DATE), minute(DATE), second(DATE), sep = ":"))
# Alle Zeilen mit Uhrzeit 23:59:00 entfernen, da letzter Tagesreport mit allen Werten NA
df17 <- df17 %>%
filter(time != "23:59:0")
df17$time <- NULL
df17$DATE <- NULL
# Temperatur und Windgeschwindigkeit umrechnen
df17 <- df17 %>%
mutate(temp = round((HourlyDryBulbTemperature-32)*5/9, 2),
wdsp = round(HourlyWindSpeed*1.609344, 2))
df17$HourlyDryBulbTemperature <- NULL
df17$HourlyWindSpeed <- NULL
# Niederschlagsmenge umbenennen, Werte mit Buchstaben korrigieren, T durch 1000 ersetzten, in mm umrechnen und die ehemaligen T-Werte in 0.01mm ändern
df17 <- df17 %>%
rename(precip = HourlyPrecipitation)
df17$precip <- str_replace_all(df17$precip, c("s" = ""))
df17$precip[df17$precip == "T"] <- 1000
df17$precip <- as.double(df17$precip)
df17 <- df17 %>%
mutate(precip = round(precip*25.4, 2))
max(df17$precip, na.rm = T)
df17$precip[df17$precip == 25400] <- 0.01
# Reihenfolge der Spalten ändern und nur eine Zeile pro Stunde behalten
df17 <- df17[,c(2,3,4,5,1)]
df17 <- df17 %>%
distinct(date, hour, .keep_all = T)
# Fehlende Werte interpolieren
df17[,3:5] <- na_interpolation(df17[,3:5], option = "spline")
apply(df17, 2, function(x) sum(is.na(x)))
# Daten speichern
saveRDS(df16, "Data/hourly_weather_2016.RDS")
saveRDS(df17, "Data/hourly_weather_2017.RDS")#Daten laden
hourly_weather_16 <- read_rds("Data/hourly_weather_2016.RDS")
hourly_weather_17 <- read_rds("Data/hourly_weather_2017.RDS")Im weiteren Verlauf dieser Arbeit wird mit den Datensätzen hourly_weather_16 und hourly_weather_2017 gearbeitet.
Daten Aggregieren
Stündliche Daten
Um die erste Forschungsfrage beantworten zu können, müssen die Daten zunächst so aggregiert werden, dass die Anzahl der Fahrten pro Stunde und Station sichtbar werden. Hier werden zwei Datensätze erstellt: Der Erste hält die Anzahl der pro Stunde losgefahrenen Radfahrer je Station fest (hourly_starts). Der Zweite hält die Anzahl der pro Stunde ankommenden Radfahrer je Station fest (hourly_stops). Die Anzahl der Fahrer wird mit der Variable user_count beschrieben. Die Variablen avg_age und avg_tripduration beschreiben jeweils das durchschnittliche Alter in Jahren, sowie die durchschnittliche Fahrtdauer in Sekunden. Die Variable weekend gibt an, ob der jeweilige Tag ein Wochenendtag (Samstag oder Sonntag) ist (1 = ja, 0 = nein). Die Variablen male_user_count, female_user_count und undefined_user_count geben die Anzahl der Fahrer nach Geschlecht gesplittet an. subscriber_count und customer_count geben an, wie viele der Fahrer ein Abo hatten bzw. einen 1-Tages oder 3-Tages-Pass. In den letzten 5 genannten Variablen werden fehlende Werte durch 0 ersetzt, da sich diese durch das Pivotieren der jeweiligen Data Frames ergeben und für die Anzahl der Fahrer stehen.
# Datensatz mit getrenntem Datum und Stunde erstellen
bike_q1 <- bike
bike_q1 <- bike_q1 %>%
mutate(start_date = as_date(starttime),
start_hour = hour(starttime),
stop_date = as_date(stoptime),
stop_hour = hour(stoptime))
# Datensatz für stündliche Starts je Station erstellen und Spalten für die Anzahl der losfahrenden Radfahrer, deren durchschnittliches Alter und die durchschnittliche Fahrtdauer einfügen, sowie Variable für Wochenende
hourly_starts <- bike_q1 %>%
group_by(start_date, start_hour, start_station_id, start_station_name, start_station_latitude,
start_station_longitude) %>%
summarise(user_count = n(),
avg_age = round(2016 - mean(birth_year)),
avg_tripduration = mean(tripduration)) %>%
mutate(weekend = ifelse(wday(start_date) == 6 | wday(start_date) == 7, 1,0))
# Datensatz erstellen der nach Stunde, Station und Geschlecht getrennt Fahrer zählt
gender <- bike_q1 %>%
group_by(start_date, start_hour, start_station_id, gender) %>%
summarise(user = n())
gender <- pivot_wider(gender, names_from = gender, values_from = user)
# Beide Datensätze verbinden und redundante Spalten löschen
hourly_starts <- cbind(hourly_starts, gender)
hourly_starts$start_date1 <- NULL
hourly_starts$start_hour1 <- NULL
hourly_starts$start_station_id1 <- NULL
# Datensatz erstellen, der nach Stunde, Station und Nutzertyp getrennt Fahrer zählt
usertype <- bike_q1 %>%
group_by(start_date, start_hour, start_station_id, usertype) %>%
summarise(user = n())
usertype <- pivot_wider(usertype, names_from = usertype, values_from = user)
# Beide Datensätze verbinden und redundante Spalten löschen
hourly_starts <- cbind(hourly_starts, usertype)
hourly_starts$start_date1 <- NULL
hourly_starts$start_hour1 <- NULL
hourly_starts$start_station_id1 <- NULL
# Spalten umbenennen
hourly_starts <- hourly_starts %>%
rename(male_user_count = "1",
female_user_count = "2",
undefined_user_count = "0",
subscriber_count = Subscriber,
customer_count = Customer)
# Fehlende Werte ersetzen
hourly_starts$male_user_count[is.na(hourly_starts$male_user_count)] <- 0
hourly_starts$female_user_count[is.na(hourly_starts$female_user_count)] <- 0
hourly_starts$undefined_user_count[is.na(hourly_starts$undefined_user_count)] <- 0
hourly_starts$subscriber_count[is.na(hourly_starts$subscriber_count)] <- 0
hourly_starts$customer_count[is.na(hourly_starts$customer_count)] <- 0
# Nicht mehr benötigte Datensätze löschen
rm(gender)
rm(usertype)
# Datensatz für stündliche Stops je Station erstellen und Spalten für die Anzahl der losfahrenden Radfahrer, deren durchschnittliches Alter und die durchschnittliche Fahrtdauer einfügen, sowie Variable für Wochenende
hourly_stops <- bike_q1 %>%
group_by(stop_date, stop_hour, end_station_id, end_station_name, end_station_latitude,
end_station_longitude) %>%
summarise(user_count = n(),
avg_age = round(2016 - mean(birth_year)),
avg_tripduration = mean(tripduration)) %>%
mutate(weekend = ifelse(wday(stop_date) == 6 | wday(stop_date) == 7, 1,0))
# Datensatz erstellen der nach Stunde, Station und Geschlecht getrennt Fahrer zählt
gender <- bike_q1 %>%
group_by(stop_date, stop_hour, end_station_id, end_station_name, end_station_latitude,
end_station_longitude,gender) %>%
summarise(user = n())
gender <- pivot_wider(gender, names_from = gender, values_from = user)
# Beide Datensätze verbinden und redundante Spalten löschen
hourly_stops <- cbind(hourly_stops, gender)
hourly_stops$stop_date1 <- NULL
hourly_stops$stop_hour1 <- NULL
hourly_stops$end_station_id1 <- NULL
hourly_stops$end_station_name1 <- NULL
hourly_stops$end_station_latitude1 <- NULL
hourly_stops$end_station_longitude1 <- NULL
# Datensatz erstellen, der nach Stunde, Station und Nutzertyp getrennt Fahrer zählt
usertype <- bike_q1 %>%
group_by(stop_date, stop_hour, end_station_id, end_station_name, end_station_latitude,
end_station_longitude, usertype) %>%
summarise(user = n())
usertype <- pivot_wider(usertype, names_from = usertype, values_from = user)
# Beide Datensätze verbinden und redundante Spalten löschen
hourly_stops <- cbind(hourly_stops, usertype)
hourly_stops$stop_date1 <- NULL
hourly_stops$stop_hour1 <- NULL
hourly_stops$end_station_id1 <- NULL
hourly_stops$end_station_name1 <- NULL
hourly_stops$end_station_latitude1 <- NULL
hourly_stops$end_station_longitude1 <- NULL
# Spalten umbenennen
hourly_stops <- hourly_stops %>%
rename(male_user_count = "1",
female_user_count = "2",
undefined_user_count = "0",
subscriber_count = Subscriber,
customer_count = Customer)
# Fehlende Werte ersetzen
hourly_stops$male_user_count[is.na(hourly_stops$male_user_count)] <- 0
hourly_stops$female_user_count[is.na(hourly_stops$female_user_count)] <- 0
hourly_stops$undefined_user_count[is.na(hourly_stops$undefined_user_count)] <- 0
hourly_stops$subscriber_count[is.na(hourly_stops$subscriber_count)] <- 0
hourly_stops$customer_count[is.na(hourly_stops$customer_count)] <- 0
# Nicht mehr benötigte Datensätze löschen
rm(gender)
rm(usertype)
# Aggregierte Datensätze speichern
saveRDS(hourly_starts, "Data/bike_q1_starts.rds")
saveRDS(hourly_stops, "Data/bike_q1_stops.rds")Tägliche Daten
Um die Daten für die zweite Frage zu aggregieren, müssen die Abfahrten und Ankünfte pro Station pro Tag zusammengefasst werden. Deswegen werde ich als erstes das Datum extrahieren, damit ich damit nach dem Tag aggregieren kann.
# Kopiere Datensatz in Question2 Datensatz und wandel Date/time Spalte in Date spalte um und füge Sie hinzu
bike_q2 <- bike
bike_q2 <- bike_q2 %>%
mutate(date_start=date(bike_q2$starttime),
date_stop=date(bike_q2$stoptime))
# Extrahiere das Gender aus dem Datensatz für die Starts
bike_q2_starts_gender <- bike_q2 %>%
group_by(date_start, start_station_name, start_station_latitude, start_station_longitude, gender) %>%
summarise( anzahl = n() )
# Extrahiere das Gender aus dem Datensatz für die Stops
bike_q2_stops_gender <- bike_q2 %>%
group_by(date_stop, end_station_name, end_station_latitude, end_station_longitude, gender) %>%
summarise( anzahl = n() )
# Tidy Gender Tabelle für starts und stops
gender_q2_starts <- pivot_wider(bike_q2_starts_gender, names_from = gender, values_from = anzahl)
gender_q2_stops <- pivot_wider(bike_q2_stops_gender, names_from = gender, values_from = anzahl)
# Aggregiere Starts und Stops nach Datum und Station
bike_q2_starts <- bike_q2 %>%
group_by(date_start, start_station_name, start_station_latitude, start_station_longitude) %>%
summarise(avr_age = round(2016-mean(birth_year)), avg_tripduration = mean(tripduration), anzahl = n() )
bike_q2_stops <- bike_q2 %>%
group_by(date_stop, end_station_name, end_station_latitude, end_station_longitude) %>%
summarise(avr_age = round(2016-mean(birth_year)), avg_tripduration = mean(tripduration), anzahl = n() )
# Füge den agregierten Tabellen das Gender hinzu
bike_q2_starts_rdy <- left_join(bike_q2_starts, gender_q2_starts, by = c("date_start", "start_station_name" ))
bike_q2_stops_rdy <- left_join(bike_q2_stops, gender_q2_stops, by = c("date_stop", "end_station_name" ))
# Nicht benötigte Spalten gelöscht und NAs in 0 umgewandelt weil wir wissen das diese 0 sind
bike_q2_starts_rdy$start_station_latitude.y <- NULL
bike_q2_starts_rdy$start_station_longitude.y <- NULL
bike_q2_starts_rdy[is.na(bike_q2_starts_rdy$`1`),"1"] <- 0
bike_q2_starts_rdy[is.na(bike_q2_starts_rdy$`2`),"2"] <- 0
bike_q2_starts_rdy[is.na(bike_q2_starts_rdy$`0`),"0"] <- 0
bike_q2_starts_rdy <- rename(bike_q2_starts_rdy, male = "1")
bike_q2_starts_rdy <- rename(bike_q2_starts_rdy, woman = "2")
bike_q2_starts_rdy <- rename(bike_q2_starts_rdy, undefined = "0")
bike_q2_stops_rdy$end_station_latitude.y <- NULL
bike_q2_stops_rdy$end_station_longitude.y <- NULL
bike_q2_stops_rdy[is.na(bike_q2_stops_rdy$`1`),"1"] <- 0
bike_q2_stops_rdy[is.na(bike_q2_stops_rdy$`2`),"2"] <- 0
bike_q2_stops_rdy[is.na(bike_q2_stops_rdy$`0`),"0"] <- 0
bike_q2_stops_rdy <- rename(bike_q2_stops_rdy, male = "1")
bike_q2_stops_rdy <- rename(bike_q2_stops_rdy, woman = "2")
bike_q2_stops_rdy <- rename(bike_q2_stops_rdy, undefined = "0")
rm(bike_q2, bike_q2_starts, bike_q2_starts_gender, bike_q2_stops, bike_q2_stops_gender, gender_q2_starts, gender_q2_stops)
saveRDS(bike_q2_starts_rdy, "Data/bike_q2_starts.rds")
saveRDS(bike_q2_stops_rdy, "Data/bike_q2_stops.rds")Visualisierung
Basis Visualisierungen
Um ein Gefühl für die Daten zu bekommen und erste Einblicke in deren Beschaffenheit zu erhalten, werden einfache Visualisierungen genutzt.
Zunächst werden die monatlichen Nutzerzahlen betrachtet.
png("Nutzerzahlen_Pro_Monat.png")
q2_starts %>%
group_by(month = month(date, label = T)) %>%
summarise(user_sum = sum(anzahl)) %>%
ggplot(aes(month, user_sum)) +
geom_line(colour = "cadetblue", size = 0.75, group = 1) +
geom_point(colour = "cadetblue", size = 1.5) +
labs(title = "Nutzeranzahl Pro Monat", y = "Nutzeranzahl",
x = "Monat")
dev.off()Nutzerzahlen Pro Monat
Diese entsprechen weitestgehend den Erwartungen. In den Wintermonaten sind die Nutzerzahlen wesentlich geringer als in den Sommermonaten. Ab Februar steigen die Nutzerzahlen bis September stetig an, mit Ausnahme eines Einbruchs im Juli. Dieser kann verschiedene Gründe haben. Zum Einen könnten viele Menschen ihren Sommerurlaub nehmen und verreisen oder die Temperaturen könnten so hoch sein, dass Radfahren als unangenehm empfunden wird.
Als nächstes sollen die monatlichen Nutzerzahlen von Männern und Frauen verglichen werden, um zu sehen, ob Unterschiede bestehen.
png("Nutzerzahlen_Pro_Monat_Nach_Geschlecht.png")
q2_starts %>%
group_by(month = month(date, label = T)) %>%
summarise(male_sum = sum(male),
female_sum = sum(woman)) %>%
ggplot() +
geom_line(aes(month, male_sum, colour = "cadetblue"), size = 0.75, group = 1) +
geom_line(aes(month, female_sum, colour = "cadetblue4"), size = 0.75, group = 1) +
scale_colour_manual(name = "Gender", values = c("cadetblue" = "cadetblue", "cadetblue4" = "cadetblue4"),
labels = c("Male", "Female")) +
labs(title = "Nutzeranzahl Pro Monat Nach Geschlecht",
y = "Nutzeranzahl", x = "Monat")
dev.off()Nutzerzahlen Pro Monat Nach Geschlecht
Auch nach der Trennung der Nutzerzahlen nach Geschlecht ist das vorherige Muster noch sichtbar. Jedoch ist klar erkennbar, dass deutlich mehr Männer das Fahrradangebot nutzen als Frauen. Hier kann weitere Forschung seitens des Unternehmens oder der Stadt New York betrieben werden, um herauszufinden, warum im Verhältnis zu den Männern so wenige Frauen das Angebot nutzen. Hieraus können dann Strategien abgeleitet werden, um mehr Frauen anzusprechen und sie als Kundinnen zu gewinnen.
In einem nächsten Schritt werden die durchschnittlichen Nutzerzahlen, sowie die durchschnittliche Fahrtdauer pro Wochentag betrachtet.
png("Durchschnittliche_Nutzeranzahl_Pro_Wochentag.png")
q2_starts %>%
group_by(wday = wday(date, label = T, abbr = F, week_start = 1)) %>%
summarise(user_mean = round(sum(anzahl)/365)) %>%
ggplot(aes(wday, user_mean)) +
geom_line(colour = "cadetblue", size = 0.75, group = 1) +
geom_point(colour = "cadetblue", size = 1.5) +
labs(title = "Durschnittliche Nutzeranzahl Pro Wochentag", y = "Durchschnittliche Nutzeranzahl",
x = "Wochentag")
png("Durchschnittliche_Fahrtdauer_Pro_Wochentag.png")
q2_starts %>%
group_by(wday = wday(date, label = T, abbr = F, week_start = 1)) %>%
summarise(mean_tripduration = mean(avg_tripduration)) %>%
ggplot(aes(wday, mean_tripduration)) +
geom_line(colour = "cadetblue", size = 0.75, group = 1) +
geom_point(colour = "cadetblue", size = 1.5) +
labs(title = "Durschnittliche Fahrtdauer Pro Wochentag", y = "Durchschnittliche Fahrtdauer",
x = "Wochentag")
dev.off()
Betrachtet man zunächst die Nutzerzahlen pro Wochentag, so fällt auf, dass diese unter der Woche höher sind als am Wochenende. Dies lässt darauf schließen, dass viele Kuden das Rad nutzen, um zur Arbeit, zur Schule oder zu Universität zu fahren. Im Gegensatz hierzu steht die durchschnittliche Dauer einer Fahrt. Die Fahrtdauer ist mit 15 Minuten am Samstag am höchsten. Dies könnte darauf schließen lassen, dass die Fahrtwege, die zur Arbeit zurückgelegt werden etwas kürzer sind, als diejenigen, die zu Freizeitzwecken zurückgelegt werden. Es kann jedoch nicht abschließend geklärt werden, ob dies tatsächlich der Fall ist, da über die Kunden außer ihrem Alter, Geschlecht und Abo-Status nichts bekannt ist. Des Weiteren sind die Unterschiede in der Fahrtdauer mit einigen Minuten nicht sehr groß.
Um weitere Anhaltspunkte für die Nutzung zu erhalten, werden die Nutzerzahlen eines Tages in Stundenintervallen betrachtet.
png("Durchschnittliche_Nutzeranzahl_Pro_Stunde.png")
hourly_starts %>%
group_by(start_hour) %>%
summarise(user_mean = round(sum(user_count)/365)) %>%
ggplot(aes(start_hour, user_mean)) +
geom_line(colour = "cadetblue4", size = 0.75) +
geom_point(colour = "cadetblue4", size = 1.5) +
labs(title = "Durchschnittliche Nutzeranzahl Pro Stunde", x = "Stunde",
y = "Durchschnittliche Nutzeranzahl")
dev.off()Durchschnittliche Nutzeranzahl Pro Stunde
Auch hier kann ein Muster erkannt werden. Von 5 Uhr bis 9 Uhr morgens steigen die Nutzerzahlen und fallen danach ab, um von 15 bis 17 Uhr wieder anzusteigen. Dies spricht für die zuvor geäußerte Vermutung, dass viele Kunden die Fahrräder für den Arbeitsweg nutzen. Am Vormittag und Nachmittag sind die Nutzerzahlen am Höchsten. Dies kann den Arbeitsbeginn und das Arbeitsende markieren.
Als nächstes werden nun die Nutzerzahlen getrennt nach Abonnement (usertype) betrachtet, um festzustellen, ob Unterschiede in der täglichen Nutzung bestehen.
png("Durchschnittliche_Nutzeranzahl_Pro_Stunde_Nach_Abo.png")
hourly_starts %>%
group_by(start_hour) %>%
summarise(subscriber_mean = round(sum(subscriber_count)/365),
customer_mean = round(sum(customer_count)/365)) %>%
ggplot() +
geom_line(aes(start_hour, subscriber_mean, colour = "cadetblue"), size = 0.75) +
geom_line(aes(start_hour, customer_mean, colour = "cadetblue4"), size = 0.75) +
scale_colour_manual(name = "User Type", values = c("cadetblue" = "cadetblue", "cadetblue4" = "cadetblue4"),
labels = c("Subscriber", "Customer")) +
labs(title = "Durchschnittliche Nutzeranzahl Pro Stunde Nach Abonnement",
y = "Durchschnittliche Nutzeranzahl", x = "Stunde")
dev.off()Durchschnittliche Nutzeranzahl Pro Stunde Nach Abonnement
Es ist zu sehen, dass für die Nutzergruppe, die ein jährliches Abo hat (Subscriber), das stündliche Muster erhalten bleibt. Jeweils morgens und nachmittags zu Arbeitsbeginn und -ende nutzen die meisten Kunden das Angebot. Die Nutzergruppe Customer, die nur einen ein- oder drei-Tages-Pass nutzen ist stark unterrepräsentiert. Dies kann darauf schließen lassen, dass das Angebot eher von Personen genutzt wird, die dies regelmäßig nutzen wollen, wie z.B. für den Arbeitsweg, und eher nicht für z.B. Tagesausflüge. Ähnlich wie bei den starken Unterschieden in der Nutzung bei den Geschlechtern, kann es hier lohnenswert sein die Hintergründe zu betrachten und weiter zu erforschen, um neue Nutzergruppen zu aquirieren und das Angebot für vorhandene Nutzergruppen zu optimieren.
Um mehr über die Nutzer zu erfahren, wird nun das durchschnittliche Alter der Nutzer pro Stunde betrachtet.
png("Durchschnittliches_Alter_Der_Nutzer_Pro_Stunde.png")
hourly_starts %>%
group_by(start_hour) %>%
summarise(mean_age = mean(avg_age)) %>%
ggplot(aes(start_hour, mean_age)) +
geom_line(colour = "cadetblue", size = 0.75) +
geom_point(colour = "cadetblue", size = 1.5) +
labs(title = "Durschnittliches Alter Der Nutzer Pro Stunde", y = "Durchschnittliches Alter",
x = "Stunde")
dev.off()Durchschnittliches Alter Der Nutzer Pro Stunde
Es ist zu erkennen, dass die Nutzer am Abend und in der Nacht jünger sind, als die Nutzer am Tag. Hierfür kann es verschiedene Gründe geben. Zum Beispiel kann es sein, dass jüngere Kunden eher nachts ausgehen und für den Heimweg das Rad nehmen. Dies kann allerdings nicht abschließend bestätig werden, da zu wenig über die Nutzer bekannt ist. Es kann sich hier ebenfalls lohnen, die Altersstruktur der Nutzer genauer zu betrachten, um neue Nutzergruppen zu identifizieren. Es könnte beispielsweise sein, dass sich die Nutzer, die mit dem Rad zu Arbeit fahren, in mehr als der Nutzungsart von denjenigen unterscheiden, die nachts mit dem Rad von einer Feier nach Hause fahren. In der Konsequenz können die unterschiedlichen Nutzergruppen besser abgestimmt auf ihr Alter und den Nutzungszweck angesprochen werden.
In einer weiteren Grafik wird die durchschnittliche Fahrtdauer pro Stunde betrachtet.
png("Durchschnittliche_Fahrtdauer_Pro_Stunde.png")
hourly_starts %>%
group_by(start_hour) %>%
summarise(mean_tripduration = mean(avg_tripduration)) %>%
ggplot(aes(start_hour, mean_tripduration)) +
geom_line(colour = "cadetblue", size = 0.75) +
geom_point(colour = "cadetblue", size = 1.5) +
labs(title = "Durschnittliche Fahrtdauer Pro Stunde", y = "Durchschnittliche Fahrtdauer",
x = "Stunde")
dev.off()Durchschnittliche Fahrtdauer Pro Stunde
Die durchschnittliche Fahrtzeit pro Stunde unterscheidet sich nur leicht. Sie bewegt sich in einem Rahmen von 700 bis 900 Sekunden bzw. 11.5 und 15 Minuten. Am längsten ist die Fahrtzeit in der Zeit zwischen 2 und 3 Uhr morgens. Dies könnte darauf hindeuten, dass die Kunden nachts eher weiter entfernt von ihrem Wohnort aufhalten.
In der letzten Grafik werden die Nutzerzahlen je Durchschnittstemperatur betrachtet.
png("Nutzeranzahl_Je_Durchschnittstemperatur.png")
q2_starts %>%
full_join(weather, by = c("date" = "date")) %>%
group_by(average.temperature) %>%
summarise(user_sum = sum(anzahl)) %>%
ggplot(aes(x = average.temperature, y = user_sum)) +
geom_point(colour = "cadetblue3", size = 1.5) +
geom_smooth(colour = "cadetblue4", se = F) +
labs(title = "Nutzeranzahl je Durchschnittstemperatur", x = "Durchschnittstemperatur in Grad Celsius",
y = "Nutzeranzahl") +
ylim(0,400000)
dev.off()Nutzeranzahl Je Durchschnittstemperatur
Es ist eine klare Tendenz erkennbar. Mit steigender Temperatur steigt die Anzahl der Nutzer. Ab einer gewissen Temperatur sinken die Nutzerzahlen wieder. Hier lässt sich vermuten, dass es den Nutzern zu warm ist, um Fahrrad zu fahren.
Ein Tag in New York
Um die Nutzerströme zwischen den Stationen zu den verschiedenen Tageszeiten sichtbar zu machen, folgt eine animierte Grafik, die jeden Nutzer zeigt, der am 01.06.2016 ein Rad entliehen und wieder abgegeben hat.
# Daten für Juni 2016 laden
df6 <- read_rds("Data/Bike_Data_Jun.RDS")
df6$starttime <- mdy_hms(df6$starttime)
df6$stoptime <- mdy_hms(df6$stoptime)
# Spalte track_id hinzufügen und Daten so aufbereiten, dass jeder Start und Stopp eine eigene Zeile hat
df6 <- rowid_to_column(df6, "track_id")
geo <- df6 %>%
pivot_longer(cols = c(starttime, stoptime), names_to = "datetime", values_to = "timestamp")
geo1 <- df6 %>%
pivot_longer(cols = c(start_station_id, end_station_id), names_to = "station", values_to = "station_id")
geo2 <- df6 %>%
pivot_longer(cols = c(start_station_latitude, end_station_latitude), names_to = "station", values_to = "station_lat")
geo3 <- df6 %>%
pivot_longer(cols = c(start_station_longitude, end_station_longitude), names_to = "station", values_to = "station_long")
geo_rdy <- cbind(geo, geo1[,c("station_id")], geo2[,c("station_lat")], geo3[,c("station_long")])
geo_rdy <- as_tibble(geo_rdy)
geo_final <- geo_rdy %>% dplyr::select(track_id, timestamp, station_id, station_lat, station_long)
# Nicht mehr genutzte Datensätze entfernen
rm(geo)
rm(geo1)
rm(geo2)
rm(geo3)
rm(geo_rdy)
# Nur komplette Fälle vom 01.06.2016 behalten, die unterschiedliche Start- und Endstationen haben
geo_final <- geo_final %>%
filter(date(timestamp) == "2016-06-01")
geo_unique <- unique(c(geo_final[,1], geo_final[,3], geo_final[,4], geo_final[,5]))
geo_unique <- as.data.frame(geo_unique)
colnames(geo_unique) <- c("track_id", "station_id", "station_lat", "station_long")
geo_unique <- geo_unique %>% distinct()
geo_unique <- geo_unique[geo_unique$track_id %in% names(which(table(geo_unique$track_id) == 2)), ]
geo_final <- inner_join(geo_unique, geo_final, by = c("track_id", "station_id", "station_lat", "station_long"))
# Karte importieren, auf die projiziert werden soll
map <- geojson_read("Data/new-york-city-boroughs.geojson", what = "sp")
# Objekt der Klasse move erstellen
geo_move <- df2move(geo_final, proj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",
track_id = "track_id", x = "station_long", y = "station_lat", time = "timestamp")
# Gemeinsame Skala für Bewegungen erstellen
geo_aligned <- align_move(geo_move, unit = "hours")
# Frames erstellen, die Bewegung zeigen
geo_frames <- frames_spatial(geo_aligned, map_service = "osm", map_type = "watercolor", alpha = 0.5) %>%
add_labels(x = "Longitude", y = "Latitude", title = "Movement Of Cyclists",
subtitle = "First Week Of June 2016") %>%
add_northarrow() %>%
add_scalebar() %>%
add_timestamps(geo_aligned, type = "label") %>%
add_progress()
# Frames animieren und als .gif speichern
animate_frames(geo_frames, out_file = "movement_of_cyclists_6_2016.gif")Die hier erstellte Grafik wird nicht weiter genutzt und betrachtet. Dies hat verschiedene Gründe. Zum Einen ist die Visualisierung aufgrund der Menge an Datenpunkten, die sich auf einem relativ kleinen Ausschnitt der Karte befinden zu unübersichtlich. Zum Anderen ist auch kein klares Fahrtenmuster erkennbar. Dies zum Teil aufgrund der Unübersichtlichkeit, aber auch, da fast alle Stationen im Stadtteil Manhattan liegen und die Fahrten recht gleichmäßig auf alle Stationen verteilt sind.
Heat Map / Choropleth Map / Interactive Map / Animierte Map
Als erstes laden wir einige Plugins zum visualisieren von Maps und Geodaten.
# load packages
library(tigris)
library(leaflet)
library(sp)
library(maptools)
library(broom)
library(httr)
library(rgdal)Nun kreiere ich einen Stations Dataframe mit den Gesamtnutzerzahlen für das gesamte Jahr.
# create stations data frame
stations <- q2_starts %>%
group_by(start_station_name, start_station_latitude.x, start_station_longitude.x) %>%
summarise( count = n(), totalusers = sum(anzahl))
stations$count <- NULLMit der GET funktion hole ich mir die verschiedenen Neighboorhoods als GeoJSON von einem Link und speicherere diese als GEOJSON.
# get the neighborhoods
r <- GET('http://data.beta.nyc//dataset/0ff93d2d-90ba-457c-9f7e-39e47bf2ac5f/resource/35dd04fb-81b3-479b-a074-a27a37888ce7/download/d085e2f8d0b54d4590b1e7d1f35594c1pediacitiesnycneighborhoods.geojson')
nyc_neighborhoods <- readOGR(content(r,'text'), 'OGRGeoJSON', verbose = F)
saveRDS(nyc_neighborhoods, "Data/GEOJSON.rds")GeoJSON laden, umformen für GGPLOT und einmal ausgeben.
nyc_neighborhoods <- readRDS("Data/GEOJSON.rds")
# tidy neighborhoods von spatial polygon data frame in ein plain data frame um es mit ggplot zu nutzen
nyc_neighborhoods_df <- tidy(nyc_neighborhoods)
ggplot() +
geom_polygon(data=nyc_neighborhoods_df, aes(x=long, y=lat, group=group), color="blue", fill=NA)Ausgabe der Spatial Data mit Leaflet.
# benutze leaflet mit neighbourhoods
leaflet(nyc_neighborhoods) %>%
addTiles() %>%
addPolygons(popup = ~neighborhood) %>%
addProviderTiles("CartoDB.Positron")Einmal Stations als points definieren und in einen spatial polygon dataframe umformen.
# lade stations als points
points <- stations %>%
dplyr::select(start_station_longitude.x, start_station_latitude.x, totalusers)
points$start_station_name <- as.factor(points$start_station_name)
points <- as.data.frame(points)
# formatiere points in ein spatial polygon data frame
points_spdf <- points
coordinates(points_spdf) <- ~start_station_longitude.x + start_station_latitude.x
proj4string(points_spdf) <- proj4string(nyc_neighborhoods)
matches <- over(points_spdf, nyc_neighborhoods)
points <- cbind(points, matches)Einmal die Karte angucken mit leaflet und eingetragenene Stationen und Neighborhoods mit den Nutzerzahlen.
# group by neighborhood
points_by_neighborhood <- points %>%
group_by(neighborhood) %>%
summarize(userperneighborhood=sum(totalusers))
map_data <- geo_join(nyc_neighborhoods, points_by_neighborhood, "neighborhood", "neighborhood")
pal <- colorNumeric(palette = "RdBu",
domain = range(map_data@data$userperneighborhood, na.rm=T))
# Karte angucken mit Stationen und Neighborhoods
leaflet(map_data) %>%
addTiles() %>%
addPolygons(fillColor = ~pal(userperneighborhood), popup = ~neighborhood) %>%
addMarkers(~start_station_longitude.x, ~start_station_latitude.x, popup = ~start_station_name, data = points) %>%
addProviderTiles("CartoDB.Positron") %>%
setView(-73.98, 40.75, zoom = 13) %>%
addLegend( pal=pal, values=~userperneighborhood, opacity=0.9, title = "Totalusers", position = "bottomright" )Die Karte ohne die Stationen mit nur den Neighborhoods und der Einfärbung der Neighborhoods nach Gesamtnutzerzahlen erstellen.
# map ohne marker
leaflet(map_data) %>%
addTiles() %>%
addPolygons(fillColor = ~pal(userperneighborhood), popup = ~neighborhood) %>%
addProviderTiles("CartoDB.Positron") %>%
setView(-73.98, 40.75, zoom = 13) %>%
addLegend( pal=pal, values=~userperneighborhood, opacity=0.9, title = "Totalusers", position = "bottomright" )Eine Karte mit Stationen ohne Neighborhoods aber mit Einfärbung der Stationen nach Nutzerzahlen erstellen.
# Eigene color Palette mit eigenen Reichweiten
mybins <- seq(1000, 50000, by=5000)
mypalette <- colorBin( palette="YlOrBr", domain=stations$totalusers, na.color="transparent", bins=mybins)
# Text für Tooltip preperieren
mytext <- paste(
"Name: ", stations$start_station_name, "<br/>",
"Totalusers: ", stations$totalusers, sep="") %>%
lapply(htmltools::HTML)
# Map kreieren
m <- leaflet(stations) %>%
addTiles() %>%
setView( lat=40, lng=-74 , zoom=4) %>%
addProviderTiles("Esri.WorldImagery") %>%
addCircleMarkers(~start_station_longitude.x, ~start_station_latitude.x, fillColor = ~mypalette(totalusers), fillOpacity = 0.7, color="white", radius=8, stroke=FALSE,
label = mytext,
labelOptions = labelOptions( style = list("font-weight" = "normal", padding = "3px 8px"), textsize = "13px",
direction = "auto")) %>%
addLegend( pal=mypalette, values=~totalusers, opacity=0.9, title = "Totalusers", position = "bottomright" )
# Map laden
m Ich versuche nun die Nutzerzahlen übers Jahr in den verschiedenen Boroughs zu animieren.
library(leaflet.minicharts)
basemap <- leaflet(map_data) %>%
addTiles() %>%
addPolygons(popup = ~neighborhood) %>%
addProviderTiles("CartoDB.Positron")
basemap %>%
addMinicharts(
q2_starts$start_station_longitude.x, q2_starts$start_station_latitude.x,
type = "auto",
chartdata = q2_starts$anzahl,
time = q2_starts$date_start,
colorPalette = colors,
transitionTime = 0.1,
) %>%
setView(-73.98, 40.75, zoom = 13)Bei der Animation animiert er nur den ersten Tag und danach keinen weiteren Tag. Das Problem wurde in der nächsten Animation durch andere Aggregation gelöst.
Neuer versuch mit Daten, die nach den Neighborhoods aggregiert wurden. Dazu aggregiere ich erstmal einen anderen Datensatz.
newTest <- bike
newTest <- newTest %>%
mutate(start_date = as_date(starttime))
newTest <- newTest %>%
left_join(points, by = c("start_station_name"="start_station_name",
"start_station_latitude"="start_station_latitude.x",
"start_station_longitude"="start_station_longitude.x"))
newTest <- newTest %>%
mutate(stop_date=date(newTest$stoptime))
newTest$totalusers <- NULL
newTest$X.id <- NULL
newTest$usertype <- NULL
newTest$start_hour <- NULL
newTest$starttime <- NULL
newTest$stoptime <- NULL
newTest$boroughCode <- NULL
newTest$borough <- NULL
newTest$neighborhood[newTest$start_station_id==160 & is.na(newTest$neighborhood)] <- "Murray Hill"
newTest$neighborhood <- fct_explicit_na(newTest$neighborhood)
newTest_agg1 <- newTest %>%
group_by(neighborhood) %>%
summarise(long = median(start_station_longitude), lat = median(start_station_latitude))
newTest <- left_join(newTest, newTest_agg1, by = "neighborhood")
newTest_agg <- newTest %>%
group_by(neighborhood, start_date, long, lat) %>%
summarise(users = n())
saveRDS(newTest_agg, "Data/newTest_agg.rds")Diese Daten speichere ich in einer RDS und lade sie wieder damit ich den Codechunk nicht immer ausführen muss.
Load Data.
Nun folgt eine Animation der Nutzerzahlen der einzelnen Neighborhoods mit dem neuen Datensatz.
Aufbereitung Für Die Analyse
Stündliche Daten
Um die erste Forschungsfrage beantworten zu können, wird der Datensatz hourly_starts noch weiter angepasst, sowie ein weiterer Datensatz erstellt, der die stündlichen Starts nach den Bezirken (neighborhoods) gruppiert und mehrere Stationen zusammenfasst.
Zunächst wird ein Datensatz erstellt, der nach Bezirken gruppiert ist. Hierfür wird der bike Datensatz mit dem Datensatz gejoined, der die Stationen und Bezirke enthält. Aufgrund von Überschneidungen in den Längen- und Breitengraden ist für die Station Nr. 160 in einigen Stunden kein Bezirk vorhanden. Dieser ist aber bekannt und wird eingesetzt. Hiernach gibt es noch 3 fehlende Bezirkswerte. Diese gehören zu einer Station ohne Koordinaten und werden mit dem Bezirk “missing” ersetzt um explizite fehlende Werte zu erhalten. Der vorhandene Datensatz wird nach Bezirken gruppiert und so aufbereitet wie im Abschnitt Preprocessing. Hiernach wird die Variable weekday für den Wochentag erstellt, sowie die Variable Holiday, die angibt, ob es sich um einen Feiertag handelt. Es folgt ein Join mit dem Datensatz hourly_weather_16, der stündliche Wetterdaten enthält. Da zwischen dem 23.01. und 26.01.2016 keine Fahrten stattfinden, werden die entsprechenden Zeilen mit den Wetterdaten gelöscht. Es fehlen auch Wetterdaten für einige Stunden, die mit Hilfe der Spline-Interpolation sinnvoll geschätzt und ersetzt werden können.
Die zuvor eingefügte Variable Wochenende wird wieder aus dem Datensatz entfernt, da das Muster, dass am Wochenende wesentlich weniger Nutzer die Räder fahren bereits über die Variable Weekday abgedeckt wird. Dadurch wird dem Problem der Kollinearität vorgebeugt. Zu guter Letzt wird der Datentyp der Spalte neighborhood von factor zu character konvertiert, da so nicht genutzte Bezirke, die als level des factors vorhanden waren, nicht länger gespeichert werden.
# Bike Daten in einen extra Data Frame laden und Datum und Stunde trennen
df <- bike
df <- df %>%
mutate(start_date = as_date(starttime),
start_hour = hour(starttime))
# Bezirke für Stationen hinzufügen
points1 <- points %>%
select(start_station_name, start_station_latitude.x, start_station_longitude.x, neighborhood)
df <- df %>%
left_join(points1, by = c("start_station_name"="start_station_name",
"start_station_latitude"="start_station_latitude.x",
"start_station_longitude"="start_station_longitude.x"))
# Fehlende Werte überprüfen
apply(df, 2, function(x) sum(is.na(x)))
# Überprüfen wo der Bezirk fehlt
df[is.na(df$neighborhood),]
# Teilweise fehlt der Bezirk für Station 160, manuell ersetzen, da bekannt
df$neighborhood[df$start_station_id==160 & is.na(df$neighborhood)] <- "Murray Hill"
# Für Station 3240 kein Bezirk bekannt, fehlende Werte explizit machen
df$neighborhood <- fct_explicit_na(df$neighborhood)
# Nach Bezirk gruppieren und Nutzerzahlen aufsummieren
hourly_starts_nh <- df %>%
group_by(start_date, start_hour, neighborhood) %>%
summarise(user_count = n()) %>%
mutate(weekend = ifelse(wday(start_date) == 6 | wday(start_date) == 7, 1,0))
# Variable für Wochentag einfügen
hourly_starts_nh <- hourly_starts_nh %>%
mutate(weekday = wday(start_date, abbr = F, label = T, week_start = 1))
hourly_starts_nh$weekday <- factor(hourly_starts_nh$weekday, ordered = FALSE)
# Variable für Feiertag einfügen
hourly_starts_nh$holiday <- ifelse(isBizday(as.timeDate(hourly_starts_nh$start_date), holidayNYSE(2016))==FALSE,1,0)
# Variable für den Monat einfügen
hourly_starts_nh$month <- month(hourly_starts_nh$start_date, abbr = F, label = T)
hourly_starts_nh$month <- factor(hourly_starts_nh$month, ordered = FALSE)
# Variable für die Woche einfügen
hourly_starts_nh$week <- week(hourly_starts_nh$start_date)
# Wetterdaten hinzufügen
hourly_starts_nh <- hourly_starts_nh %>%
full_join(hourly_weather_16, by = c("start_date" = "date", "start_hour" = "hour"))
# Fehlende Werte überprüfen
apply(hourly_starts_nh, 2, function(x) sum(is.na(x)))
hourly_starts_nh[is.na(hourly_starts_nh$neighborhood),] #vom 23.01.-26.01. keine fahrten
# Nicht genutzte Wetterdaten entfernen
hourly_starts_nh <- hourly_starts_nh[complete.cases(hourly_starts_nh[,3:7]),]
# Fehlende Werte überprüfen und interpolieren
apply(hourly_starts_nh, 2, function(x) sum(is.na(x)))
hourly_starts_nh[,10:12] <- na_interpolation(hourly_starts_nh[,10:12], option = "spline")
# Variable Weekend entfernen, da Wochentage für Betrachtung ausreichen
hourly_starts_nh$weekend <- NULL
# Spalte Neighborhood in character umwandeln, um nicht genutzte level zu eliminieren
hourly_starts_nh$neighborhood <- as.character(hourly_starts_nh$neighborhood)
# Daten speichern
saveRDS(hourly_starts_nh, "Data/q1_starts_nh.RDS")Um stationsgenaue Vorhersagen zu treffen, wird der Datenssatz hourly_starts ebenfalls weiter aufbereitet. Wie zuvor werden die Variablen Wochentag und Feiertag eingefügt, sowie die Wetterdaten hinzugefügt.
# Variable für Wochentag einfügen
hourly_starts <- hourly_starts %>%
mutate(weekday = wday(start_date, abbr = F, label = T, week_start = 1))
hourly_starts$weekday <- factor(hourly_starts$weekday, ordered = FALSE)
# Variable für Feiertag einfügen
hourly_starts$holiday <- ifelse(isBizday(as.timeDate(hourly_starts$start_date), holidayNYSE(2016))==FALSE,1,0)
# Wetterdaten hinzufügen
hourly_starts <- hourly_starts %>%
full_join(hourly_weather_16, by = c("start_date" = "date", "start_hour" = "hour"))
# Fehlende Werte überprüfen
apply(hourly_starts_nh, 2, function(x) sum(is.na(x)))
hourly_starts[is.na(hourly_starts$start_station_name),] #vom 23.01.-26.01. keine fahrten
# Nicht mehr benötigte Variablen entfernen
hourly_starts$start_station_id <- NULL
hourly_starts$start_station_latitude <- NULL
hourly_starts$start_station_longitude <- NULL
hourly_starts$avg_age <- NULL
hourly_starts$avg_tripduration <- NULL
hourly_starts$male_user_count <- NULL
hourly_starts$female_user_count <- NULL
hourly_starts$undefined_user_count <- NULL
hourly_starts$subscriber_count <- NULL
hourly_starts$customer_count <- NULL
hourly_starts$weekend <- NULL
# Nicht genutzte Wetterdaten entfernen
hourly_starts <- hourly_starts[complete.cases(hourly_starts[,3:6]),]
# Fehlende Werte überprüfen und interpolieren
apply(hourly_starts, 2, function(x) sum(is.na(x)))
hourly_starts[,7:9] <- na_interpolation(hourly_starts[,7:9], option = "spline")Dieser Datensatz wird in den weiteren Analysen nicht betrachtet. Alle weiteren Berechnungen werden mit dem Datensatz hourly_starts_nh durchgeführt, der die Fahrerzahlen auf Neighborhood-Ebene aggregiert. Dies hat den Grund, dass der obige Datensatz die Fahrerzahlen auf Stationsebene aggregiert und nach der Dummyfizierung so groß ist, dass die vorhandene Rechenleistung für die notwenigen Berechnungen nicht ausreicht.
Auf dem Datensatz hourly_starts_nh basierend, wird ein Datensatz erstellt, der nur die stündlichen Daten des Bezirks Chelsea enthält. Da der Datensatz hourly_starts_nh sehr groß ist und während der weiteren Aufbereitung immer größer wird, sodass notwendige Berechnungen sehr lange dauern oder nicht durchführbar sind, wird ein Bezirk ausgewählt für den die Berechnungen exemplarisch durchgeführt werden. Der Bezirk Chelsea wird gewählt, da er die meisten Fahrten und beinah lückenlos für jede Stunde Daten aufweist. Die Variablen Monat und Woche werden aus diesem Datensatz entfernt, da sie im späteren Verlauf der Analyse das Modell für die Regression, sowie für das LSTM Netz stark aufblähen. Um dennoch die saisonale Komponente in den Daten zu behalten, wird eine Variable für das Quartal eingefügt. Für die spätere Visualisierung der Ergebnisse wird eine Spalte erstellt, die Datum und Stunde zusammenfasst.
Die Modellierung des stündlichen LSTM Netzes wird später ausschließlich anhand dieses Datensatzes erfolgen.
# Trainingsdaten nur für Bezirk Chelsea
hourly_chelsea_train <- hourly_starts_nh %>%
dplyr::select(-month, -week) %>%
filter(neighborhood == "Chelsea") %>%
mutate(quarter = quarter(start_date),
start_datetime = make_datetime(year = year(start_date), month = month(start_date),
day = day(start_date), hour = start_hour, min = 0, sec = 0))
hourly_chelsea_train <- hourly_chelsea_train[, c(11,1:2,4,3,10,5:9)]Tägliche Daten
Als erstes den Datensatz mit dem Wetterdatensatz joinen um einen vollständigen Datensatz zu erhalten. Dazu benennen wir die Datumsspalte in beiden Datenframes gleich und joinen sie mit einem Leftjoin.
Aggregation eines Datensatzes mit Neighborhoods. Wir joinen stations mit q2_stats_w und ersetzen den fehlenden Bezirk.
points1 <- points %>%
dplyr::select(start_station_name, start_station_latitude.x, start_station_longitude.x, neighborhood)
q2_starts_w <- left_join(q2_starts_w, points1, by = c("start_station_name"="start_station_name",
"start_station_latitude.x"="start_station_latitude.x",
"start_station_longitude.x"="start_station_longitude.x"))
#Teilweise fehlt der Bezirk für Station 160, manuell ersetzen, da bekannt
q2_starts_w$neighborhood[q2_starts_w$start_station_name=="E 37 St & Lexington Ave" & is.na(q2_starts_w$neighborhood)] <- "Murray Hill"
#Für Station 3240 kein Bezirk bekannt, fehlende Werte explizit machen
q2_starts_w$neighborhood <- fct_explicit_na(q2_starts_w$neighborhood)Als nächstes füge ich die Wochentag Variable ein und erstelle ein Holiday Feature.
#Variable für Wochentag einfügen
q2_starts_w <- q2_starts_w %>%
mutate(weekday = wday(date_start, abbr = F, label = T, week_start = 1))
q2_starts_w$weekday <- factor(q2_starts_w$weekday, ordered = FALSE)
#Variable für Feiertag einfügen
q2_starts_w$holiday <- ifelse(isBizday(as.timeDate(q2_starts_w$date_start), holidayNYSE(2016))==FALSE,1,0)Prüfen auf fehlende Werte.
Nicht brauchbare Spalten entfernen, Wochentag und Feiertag hinzufügen, Neighborhoodsspalte in character umwandeln, um Problemen vorzubeugen. Environment aufräumen und Monat und Wochentag als nicht geordnete Faktoren hinzufügen.
q2_starts_w$male <- NULL
q2_starts_w$woman <- NULL
q2_starts_w$undefined <- NULL
q2_starts_w$start_station_latitude.x <- NULL
q2_starts_w$start_station_longitude.x <- NULL
q2_starts_w2 <- q2_starts_w
q2_starts_w$neighborhood <- NULL
q2_starts_w2_agg <- q2_starts_w2 %>%
group_by(date_start, neighborhood) %>%
summarise( anzahl = sum(anzahl), avg_tripduration = mean(avg_tripduration), avg_age=mean(avr_age))
q2_starts_w2_agg <- rename(q2_starts_w2_agg,
"date" = "date_start")
q2_starts_w2_agg <- left_join(q2_starts_w2_agg, weather, by = "date")
#Variable für Wochentag einfügen
q2_starts_w2_agg <- q2_starts_w2_agg %>%
mutate(weekday = wday(date, abbr = F, label = T, week_start = 1))
q2_starts_w2_agg$weekday <- factor(q2_starts_w2_agg$weekday, ordered = FALSE)
#Variable für Feiertag einfügen
q2_starts_w2_agg$holiday <- ifelse(isBizday(as.timeDate(q2_starts_w2_agg$date), holidayNYSE(2016))==FALSE,1,0)#
#Spalte Neighborhood in character umwandeln, um nicht genutzte level zu eliminieren
q2_starts_w2_agg$neighborhood <- as.character(q2_starts_w2_agg$neighborhood)
# Remove environment
rm(q2_starts_w, q2_starts)
#Monat und Woche als Variablen eingefügt
#q2_starts_w2_agg <- q2_starts_w2_agg %>%
#mutate(week = week(date))
q2_starts_w2_agg <- q2_starts_w2_agg %>%
mutate(month = month(date, abbr = F, label = T))
#Variable für Wochentag als ungeordneten Faktor speichern
q2_starts_w2_agg$weekday <- factor(q2_starts_w2_agg$weekday, ordered = FALSE)
q2_starts_w2_agg$month <- factor(q2_starts_w2_agg$month, ordered = FALSE)
# Save Data
saveRDS(q2_starts_w2_agg, "Data/q2_starts_w2_agg.RDS")Load Data.
Testdaten
Allgemeine Aggregation der Testdaten
Um die Vorhersagegenauigkeit der nachfolgenden Modelle zu überprüfen, werden Testdaten benötigt. Da alle Modelle mit den Daten des gesamten Jahres 2016 lernen, soll mit Daten aus dem Jahr 2017 getestet werden. Um alle Jahreszeiten angemessen zu betrachten, werden die Monate Februar, Mai, August und November ausgewählt. Um Speicherplatz zu sparen, werden diese im Format .RDS gespeichert.
feb2017 <- read_csv("Data/201702-citibike-tripdata.csv")
may2017 <- read_csv("Data/201705-citibike-tripdata.csv")
aug2017 <- read_csv("Data/201708-citibike-tripdata.csv")
nov2017 <- read_csv("Data/201711-citibike-tripdata.csv")
saveRDS(feb2017, "Data/feb2017.rds")
saveRDS(may2017, "Data/may2017.rds")
saveRDS(aug2017, "Data/aug2017.rds")
saveRDS(nov2017, "Data/nov2017.rds")Ebenso wie die Daten für das Jahr 2016, müssen auch die Daten für das Jahr 2017 bearbeitet werden, bevor sie in die Analysen einfliessen können. Zunächst werden Leerzeichen in den Spaltennamen durch Unterstriche ersetzt, um Probleme beim Ansprechen der Spalten zu vermeiden. Danach werden die Spaltennamen in ein gemeinsames Format überführt.
# Datensätze laden
test_feb2017 <- read_rds("Data/feb2017.rds")
test_may2017 <- read_rds("Data/may2017.rds")
test_aug2017 <- read_rds("Data/aug2017.rds")
test_nov2017 <- read_rds("Data/nov2017.rds")
# Leerzeichen in Spaltennamen ersetzten
names(test_feb2017) <- str_replace_all(names(test_feb2017), c(" " = "_"))
names(test_may2017) <- str_replace_all(names(test_may2017), c(" " = "_"))
names(test_aug2017) <- str_replace_all(names(test_aug2017), c(" " = "_"))
names(test_nov2017) <- str_replace_all(names(test_nov2017), c(" " = "_"))
# Spaltennamen an Rest des Datensatzes anpassen
test_feb2017 <- test_feb2017 %>%
rename(tripduration = Trip_Duration,
starttime = Start_Time,
stoptime = Stop_Time,
start_station_id = Start_Station_ID,
start_station_name = Start_Station_Name,
start_station_latitude = Start_Station_Latitude,
start_station_longitude = Start_Station_Longitude,
end_station_id = End_Station_ID,
end_station_name = End_Station_Name,
end_station_latitude = End_Station_Latitude,
end_station_longitude = End_Station_Longitude,
bikeid = Bike_ID,
usertype = User_Type,
birth_year = Birth_Year,
gender = Gender)
# Einen Testdatensatz erstellen
test_data <- rbind(test_feb2017, test_may2017, test_aug2017, test_nov2017)
# Nicht mehr benötigte Objekte entfernen
rm(test_feb2017, test_aug2017, test_may2017, test_nov2017)In den Testdaten aus dem Jahr 2017 befinden sich einige Stationen mehr, als noch in 2016. Aufgrunddessen wird ein Datensatz erstellt, der die Stationen für das Jahr 2017 enthält, sowie die dazugehörigen Bezirke.
# Stationen für 2017 in einem eigenen Datensatz speichern
stations_2017 <- unique(c(test_data[,5], test_data[,6], test_data[,7]))
stations_2017 <- as.data.frame(stations_2017)
stations_2017 <- stations_2017 %>% distinct()
names(stations_2017) <- c("start_station_name", "start_station_latitude", "start_station_longitude")
# Bezirke laden
nyc_neighborhoods <- readRDS("Data/GEOJSON.rds")
# Stationen für 2017 bearbeiten
stations_2017$start_station_name <- as.factor(stations_2017$start_station_name)
# stations_2017 mit neighborhoods verbinden
points_spdf <- stations_2017
coordinates(points_spdf) <- ~start_station_longitude + start_station_latitude
proj4string(points_spdf) <- proj4string(nyc_neighborhoods)
matches <- over(points_spdf, nyc_neighborhoods)
stations_2017 <- cbind(stations_2017, matches)
# Nur die Spalten behalten, die weiterhin wichtig sind
stations_2017 <- stations_2017 %>%
dplyr::select(start_station_name, start_station_latitude, start_station_longitude, neighborhood)
# Nicht mehr benötige Objekte entfernen
rm(nyc_neighborhoods, points_spdf, matches)Da zur Beantwortung der beiden Forschungsfragen unterschiedliche Anforderungen an die Testdaten gestellt werden, werden diese einmal für die tägliche und einmal für die stündliche Analyse aufbereitet.
Aufbereitung Für Stündliche Analyse
Für die stündliche Analyse werden die Daten wie folgt aufbereitet. Wie zuvor beim Trainingsdatensatz wird die Spalte starttime in die zwei Spalten Datum und Stunde aufgesplittet. Da eine Betrachtung auf Bezirkseben erfolgt, werden die Bezirke zu dem Datensatz hinzugefügt. Für zwei Stationen ist keine Zuordnung zu einem Bezirk möglich, sodass die fehlenden Werte in explizite fehlende Werte umgewandelt werden. Dies hat zur Folge, dass der Bezirk dieser Stationen “Missing” benannt wird. Da in 2017 neue Stationen auch in Bezirken errichtet wurden, die in 2016 noch über keine Stationen verfügten, werden diese “neuen” Bezirke aus den Daten entfernt. Sie können ohne Grundlage aus dem vorherigen Jahr nicht prognostiziert werden. Wie auch die Trainingsdaten werden die Testdaten nach dem Bezirk, dem Datum und der Stunde gruppiert, um die jeweiligen Nutzerzahlen aufsummieren zu können. Die Variablen, die die Ankunftsstationen betreffen, sowie die Nutzerdaten, wie z.B. Alter, Geschlecht und Fahrtdauer, werden nicht weiter betrachtet. Sie werden nicht weiter in die Analyse einfließen, da es das Ziel ist zukünftige Fahrten zu prognostizieren. Es ist dann noch nicht bekannt welches Alter die Fahrer haben werden, wie lang sie fahren oder welches Geschlecht sie haben. Es werden ebenfalls Variablen für den Wochentag und die Feiertage eingefügt, sowie Wetterdaten angefügt.
# Einen Testdatensatz erstellen
hourly_test <- test_data
# Datum und Stunde trennen
hourly_test <- hourly_test %>%
mutate(start_date = as_date(starttime),
start_hour = hour(starttime))
# Bezirke für Stationen hinzufügen
hourly_test <- hourly_test %>%
left_join(stations_2017, by = c("start_station_name"="start_station_name",
"start_station_latitude"="start_station_latitude",
"start_station_longitude"="start_station_longitude"))
# Fehlende Werte überprüfen
apply(hourly_test, 2, function(x) sum(is.na(x)))
# Überprüfen wo der Bezirk fehlt
hourly_test[is.na(hourly_test$neighborhood),]
# Der Bezirk von zwei Stationen ist nicht bekannt, fehlende Werte explizit machen
hourly_test$neighborhood <- fct_explicit_na(hourly_test$neighborhood)
# Nach Bezirk gruppieren und Nutzerzahlen aufsummieren
hourly_test <- hourly_test %>%
group_by(start_date, start_hour, neighborhood) %>%
summarise(user_count = n()) %>%
mutate(weekday = wday(start_date, abbr = F, label = T, week_start = 1),
holiday = ifelse(isBizday(as.timeDate(start_date), holidayNYSE(2017))==FALSE,1,0),
month = month(start_date, abbr = F, label = T),
week = week(start_date))
# Variablen für Wochentag und Monat als ungeordneten Faktor speichern
hourly_test$weekday <- factor(hourly_test$weekday, ordered = FALSE)
hourly_test$month <- factor(hourly_test$month, ordered = FALSE)
# Wetterdaten hinzufügen
hourly_test <- hourly_test %>%
left_join(hourly_weather_17, by = c("start_date" = "date", "start_hour" = "hour"))
# Fehlende Werte überprüfen
apply(hourly_starts_nh, 2, function(x) sum(is.na(x)))
# Spalte Neighborhood in character umwandeln, um nicht genutzte level zu eliminieren
hourly_test$neighborhood <- as.character(hourly_test$neighborhood)
# Bezirke ausfindig machen, die in 2017 vorkommen, aber nicht in 2016 und diese eliminieren
setdiff(hourly_test$neighborhood, hourly_starts_nh$neighborhood)
hourly_test <- hourly_test %>%
filter(neighborhood != "Astoria" & neighborhood != "Crown Heights" & neighborhood != "Ditmars Steinway" &
neighborhood != "Harlem" & neighborhood != "Morningside Heights" &
neighborhood != "Prospect Heights" & neighborhood != "Prospect-Lefferts Gardens")
# Daten speichern
saveRDS(hourly_test, "Data/q1_starts_test.RDS")Wie schon bei den Trainingsdaten wird auch hier ein Testdatensatz erstellt, der ausschließlich Daten aus dem Bezirk Chelsea enthält. Da die Datensätze die gleichen Attribute haben müssen, damit eine Vorhersage mit dem LSTM Netz erfolgen kann, wird auch hier der Monat und die Woche entfernt und durch das Quartal ersetzt. Für die Visualisierung der Ergebnisse wird hier ebenfalls eine Variable erstellt, die Datum und Uhrzeit zusammenfasst.
# Testdaten nur für Bezirk Chelsea
hourly_chelsea_test <- hourly_test %>%
dplyr::select(-month, -week) %>%
filter(neighborhood == "Chelsea") %>%
mutate(quarter = quarter(start_date),
start_datetime = make_datetime(year = year(start_date), month = month(start_date),
day = day(start_date), hour = start_hour, min = 0, sec = 0))
hourly_chelsea_test <- hourly_chelsea_test[, c(11,1:2,4,3,10,5:9)]Aufbereitung Für Tägliche Analyse
Um unsere Regression bzw. unser LSTM für unsere Forschungsfrage zu überprüfen, brauchen wir zusätzlich zu den Traininsdaten auch Testdaten. Diese sind aus den Monaten Februar, Mai, August und November von 2017. Im nachfolgenden Codechunk werden diese Testdaten genauso aggregiert wie oben die Traininsdaten.
# Dataframe für aggregierte Tesdaten erstellen
test_data_q2 <- test_data
# Nicht benötigte Columns entfernen
test_data_q2$gender <- NULL
test_data_q2$usertype <- NULL
#Datum trennen
test_data_q2 <- test_data_q2 %>%
mutate(start_date = as_date(starttime))
# stations 2017 umwandeln
stations_2017_test <- stations_2017
stations_2017_test$start_station_name <- as.character(stations_2017_test$start_station_name)
stations_2017_test$neighborhood <- as.character(stations_2017_test$neighborhood)
#Bezirke für Stationen hinzufügen
test_data_q2_test <- test_data_q2 %>%
left_join(stations_2017, by = c("start_station_name"="start_station_name",
"start_station_latitude"="start_station_latitude",
"start_station_longitude"="start_station_longitude"))
test_data_q2 <- test_data_q2 %>%
left_join(stations_2017_test, by = c("start_station_name"="start_station_name",
"start_station_latitude"="start_station_latitude",
"start_station_longitude"="start_station_longitude"))
# Points1 benutzen damit man nur die Stations kriegt die auch in den Trainingsdaten vorhanden sind
'points1_test <- points1
points1_test$start_station_name <- as.character(points1_test$start_station_name)
points1_test$neighborhood <- as.character(points1_test$neighborhood)
points1_test <- points1_test %>%
rename(start_station_latitude = start_station_latitude.x,
start_station_longitude = start_station_longitude.x)
#Bezirke für Stationen hinzufügen
test_data_q2 <- test_data_q2 %>%
left_join(points1_test, by = c("start_station_name"="start_station_name",
"start_station_latitude"="start_station_latitude",
"start_station_longitude"="start_station_longitude"))'
#Fehlende Werte überprüfen
#apply(test_data_q2, 2, function(x) sum(is.na(x)))
#Überprüfen wo der Bezirk fehlt
test_data_q2[is.na(test_data_q2$neighborhood),]
#Der Bezirk von zwei Stationen ist nicht bekannt, fehlende Werte explizit machen
test_data_q2$neighborhood <- fct_explicit_na(test_data_q2$neighborhood)
#Spalte Neighborhood in character umwandeln, um nicht genutzte level zu eliminieren
test_data_q2$neighborhood <- as.character(test_data_q2$neighborhood)
test_data_q2$birth_year <- as.numeric(test_data_q2$birth_year)
# Löschen der fehlenden Werte mit Birthyear
test_data_q2 <- na.omit(test_data_q2)
#Nach Bezirk gruppieren und Nutzerzahlen aufsummieren
test_data_q2 <- test_data_q2 %>%
group_by(start_date, neighborhood) %>%
summarise(anzahl = n(), avg_age = round(2016-mean(birth_year)), avg_tripduration = mean(tripduration)) %>%
mutate(weekday = wday(start_date, abbr = F, label = T, week_start = 1),
holiday = ifelse(isBizday(as.timeDate(start_date), holidayNYSE(2017))==FALSE,1,0),)
#Fehlende Werte überprüfen
apply(test_data_q2, 2, function(x) sum(is.na(x)))
#Wetter hinzufügen
test_data_q2 <- left_join(test_data_q2, weather_2017, by = c("start_date" = "date"))
# Columns wie bei den Trainingsdaten anordnen
test_data_q2 <- test_data_q2[, c(1,2,3,5,4,8,9,10,11,12,13,6,7)]
rm(stations_2017_test, test_data_q2_test, points, points1_test)
#Bezirke ausfindig machen, die in 2017 vorkommen, aber nicht in 2016 und diese eliminieren
setdiff(test_data_q2$neighborhood, q2_starts_w2_agg$neighborhood)
test_data_q2 <- test_data_q2 %>%
filter(neighborhood != "Astoria" & neighborhood != "Crown Heights" & neighborhood != "Ditmars Steinway" &
neighborhood != "Harlem" & neighborhood != "Morningside Heights" &
neighborhood != "Prospect Heights" & neighborhood != "Prospect-Lefferts Gardens")
#Monat und Woche als Variablen eingefügt
#test_data_q2 <- test_data_q2 %>%
#mutate(week = week(start_date))
test_data_q2 <- test_data_q2 %>%
mutate(month = month(start_date, abbr = F, label = T))
#Variable für Wochentag als ungeordneten Faktor speichern
test_data_q2$weekday <- factor(test_data_q2$weekday, ordered = FALSE)
test_data_q2$month <- factor(test_data_q2$month, ordered = FALSE)
#Als .RDS speichern
saveRDS(test_data_q2, "Data/test_data_q2.RDS")Aufbereitung Für Tägliche Analyse Für Ein Neighborhood
Als Erstes werden die Daten neu aggregiert. Dazu filtere ich aus den schon aggregierten Daten des Neighborhoods Chelsea raus, um uns auf diese zu konzentrieren. Danach lösche ich die Spalte Neighborhoods, um die Zahl der Features zu minimieren.
q2_starts_w2_agg <- q2_starts_w2_agg %>%
rename(
start_date = date
)
daily_neighborhood_train <- q2_starts_w2_agg %>%
filter(neighborhood == "Chelsea")
daily_neighborhood_train$neighborhood <- NULL
daily_neighborhood_test <- test_data_q2 %>%
filter(neighborhood == "Chelsea")
daily_neighborhood_test$neighborhood <- NULL
avg_anzahl <- daily_neighborhood_train %>%
dplyr::select(anzahl)
avg_anzahl$date <- NULL
avg_anzahl <- as.data.frame(avg_anzahl)
avg_anzahl %>%
summarise(avg_anzahl = mean(anzahl))## avg_anzahl
## 1 3880.058
Lineare Regression
Um eine Vergleichsmöglichkeit für die Vorhersage der KNNs zu schaffen, wird eine lineare Regression als Benchmark durchgeführt. Es werden für die beiden Forschungsfragen lineare Regressionsmodelle erstellt, die mit den Testdaten eine Vorhersage für die Nutzeranzahl treffen. Der Fehler, der hierbei gemacht wird, wird in Form des Root Mean Squared Error (RMSE) gemessen und mit dem des KNN verglichen. Dieser gibt in diesem Fall die durchschnittliche Anzahl der Nutzer an, um die das Modell daneben liegt.
Stündliche Regression
Für die stündlichen Nutzerdaten werden drei Regressionen durchgeführt. Dies ist der Fall, da sich die Granularität der Datensätze voneinander unterscheidet. Für einen Datensatz sollen die Nutzerzahlen je Bezirk bestimmt werden, für den anderen Datensatz die Nutzerzahlen je Station und für den Dritten die Nutzerzahlen nur für einen Bezirk. Die dritte Regression wird als Benchmark für das später erstellte künstliche neuronale Netz genutzt.
In einem ersten Schritt wird eine lineare Regression auf Bezirksebene durchgeführt. Die Variable start_date wird hierbei nicht berücksichtigt, da sich mit der Variable week eine Kollinearität ergibt. Diese wird weiterhin berücksichtigt.
# Regressionsmodell mit Bezirken
model_q1_nh <- as.formula(user_count ~ . -start_date)
lin_reg_q1_nh <- lm(model_q1_nh, hourly_starts_nh)
anova(lin_reg_q1_nh)## Analysis of Variance Table
##
## Response: user_count
## Df Sum Sq Mean Sq F value Pr(>F)
## start_hour 1 23765489 23765489 8835.236 < 2.2e-16 ***
## neighborhood 48 466442661 9717555 3612.671 < 2.2e-16 ***
## weekday 6 15949222 2658204 988.234 < 2.2e-16 ***
## holiday 1 3656889 3656889 1359.512 < 2.2e-16 ***
## month 11 33814466 3074042 1142.829 < 2.2e-16 ***
## week 1 111382 111382 41.408 1.237e-10 ***
## temp 1 20385993 20385993 7578.849 < 2.2e-16 ***
## wdsp 1 83747 83747 31.135 2.409e-08 ***
## precip 1 2854353 2854353 1061.156 < 2.2e-16 ***
## Residuals 316527 851411269 2690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## lm(formula = model_q1_nh, data = hourly_starts_nh)
##
## Residuals:
## Min 1Q Median 3Q Max
## -178.59 -20.50 -3.77 13.40 1211.36
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -32.377044 29.950792 -1.081 0.279694
## start_hour 1.596521 0.014151 112.824 < 2e-16 ***
## neighborhoodBattery Park City 30.574609 29.951256 1.021 0.307344
## neighborhoodBedford-Stuyvesant 14.798440 29.951009 0.494 0.621244
## neighborhoodBoerum Hill 1.074560 29.952204 0.036 0.971381
## neighborhoodBrooklyn Heights 10.532346 29.951495 0.352 0.725103
## neighborhoodCarroll Gardens 3.680145 29.961980 0.123 0.902244
## neighborhoodCentral Park 21.123775 29.951437 0.705 0.480644
## neighborhoodChelsea 161.927113 29.950651 5.406 6.43e-08 ***
## neighborhoodChinatown 20.242432 29.950969 0.676 0.499135
## neighborhoodCivic Center 9.095917 29.951579 0.304 0.761366
## neighborhoodClinton Hill 10.964039 29.951166 0.366 0.714318
## neighborhoodCobble Hill -4.640456 29.964079 -0.155 0.876926
## neighborhoodColumbia St -5.104182 29.955208 -0.170 0.864701
## neighborhoodDowntown Brooklyn 7.711770 29.951435 0.257 0.796812
## neighborhoodDUMBO -0.755089 29.952864 -0.025 0.979888
## neighborhoodEast Harlem 2.876974 29.959715 0.096 0.923498
## neighborhoodEast Village 98.643901 29.950655 3.294 0.000989 ***
## neighborhoodFinancial District 43.164296 29.950903 1.441 0.149538
## neighborhoodFlatiron District 29.692186 29.951166 0.991 0.321514
## neighborhoodFort Greene 19.213501 29.951149 0.641 0.521202
## neighborhoodGovernors Island -9.084588 30.005810 -0.303 0.762072
## neighborhoodGowanus -0.832403 29.962571 -0.028 0.977836
## neighborhoodGramercy 35.363007 29.951008 1.181 0.237725
## neighborhoodGreenpoint 20.470390 29.951050 0.683 0.494316
## neighborhoodGreenwich Village 54.927881 29.950799 1.834 0.066664 .
## neighborhoodHell's Kitchen 79.421283 29.950710 2.652 0.008008 **
## neighborhoodKips Bay 40.798231 29.950872 1.362 0.173145
## neighborhoodLong Island City 7.735734 29.951459 0.258 0.796194
## neighborhoodLower East Side 53.245860 29.950665 1.778 0.075440 .
## neighborhoodMidtown 161.697625 29.950678 5.399 6.71e-08 ***
## neighborhoodMurray Hill 38.503191 29.951046 1.286 0.198605
## neighborhoodNavy Yard -8.169031 29.956632 -0.273 0.785088
## neighborhoodNoHo 4.216448 29.951992 0.141 0.888049
## neighborhoodNolita 7.350898 29.951643 0.245 0.806127
## neighborhoodPark Slope 4.490135 29.952755 0.150 0.880838
## neighborhoodProspect Park -1.590645 29.965041 -0.053 0.957666
## neighborhoodRed Hook -5.768641 29.966483 -0.193 0.847348
## neighborhoodSoHo 45.179488 29.950870 1.508 0.131440
## neighborhoodSouth Slope -5.713108 29.965118 -0.191 0.848793
## neighborhoodStuyvesant Town 25.498741 29.950964 0.851 0.394576
## neighborhoodSunset Park -17.889665 30.026918 -0.596 0.551318
## neighborhoodTheater District 25.925479 29.951063 0.866 0.386713
## neighborhoodTribeca 40.076475 29.950927 1.338 0.180874
## neighborhoodTwo Bridges 1.320849 29.952261 0.044 0.964826
## neighborhoodUpper East Side 65.114306 29.950815 2.174 0.029703 *
## neighborhoodUpper West Side 67.902385 29.950825 2.267 0.023383 *
## neighborhoodVinegar Hill -0.843232 29.952682 -0.028 0.977541
## neighborhoodWest Village 64.215558 29.950737 2.144 0.032030 *
## neighborhoodWilliamsburg 47.658821 29.950780 1.591 0.111557
## weekdayDienstag 1.694371 0.356513 4.753 2.01e-06 ***
## weekdayMittwoch 2.391670 0.355035 6.736 1.63e-11 ***
## weekdayDonnerstag 1.198441 0.353616 3.389 0.000701 ***
## weekdayFreitag -1.246903 0.350134 -3.561 0.000369 ***
## weekdaySamstag 4.543729 0.666768 6.815 9.47e-12 ***
## weekdaySonntag 2.687937 0.667564 4.026 5.66e-05 ***
## holiday -17.852841 0.639749 -27.906 < 2e-16 ***
## monthFebruar 1.162120 0.603868 1.924 0.054298 .
## monthMärz 4.676641 0.818224 5.716 1.09e-08 ***
## monthApril 8.238749 1.094478 7.528 5.18e-14 ***
## monthMai 9.243367 1.402555 6.590 4.39e-11 ***
## monthJuni 12.874937 1.704095 7.555 4.19e-14 ***
## monthJuli 7.819132 2.015035 3.880 0.000104 ***
## monthAugust 14.597768 2.328460 6.269 3.63e-10 ***
## monthSeptember 29.861510 2.624208 11.379 < 2e-16 ***
## monthOktober 43.355630 2.931681 14.789 < 2e-16 ***
## monthNovember 46.226275 3.246522 14.239 < 2e-16 ***
## monthDezember 49.376731 3.568117 13.838 < 2e-16 ***
## week -0.907800 0.073679 -12.321 < 2e-16 ***
## temp 1.677223 0.019238 87.181 < 2e-16 ***
## wdsp 0.039288 0.006483 6.060 1.36e-09 ***
## precip -7.322706 0.224793 -32.575 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51.86 on 316527 degrees of freedom
## Multiple R-squared: 0.3998, Adjusted R-squared: 0.3996
## F-statistic: 2969 on 71 and 316527 DF, p-value: < 2.2e-16
linreg_q1_nh_train_mse <- mean((lin_reg_q1_nh$residuals)^2)
linreg_q1_nh_train_rmse <- sqrt(mean((lin_reg_q1_nh$residuals)^2))
linreg_q1_nh_train_rmse## [1] 51.8579
Die Anova deutet darauf hin, dass der Einfluss aller Variablen auf die Anzahl an Radfahrern statistisch signifikant ist, da der p-Wert äußerst gering ist. Die Funktion lm() hat zur Durchführung der Regression Dummyvariablen für den Bezirk, den Wochentag, sowie den Monat erstellt. Sieht man sich die genaue Zusammenfassung des Modells an, wird sichtbar, dass der Einfluss vieler Bezirke statistisch nicht signifikant, d.h. zufällig ist. Sie werden dennoch nicht aus dem Modell entfernt, da sich sonst ein unvollständiges Bild der Situation ergeben würde, denn in diesen Bezirken gibt es auch Leihstationen, die genutzt werden. Alle anderen Variablen, mit Ausnahme des Monats Februar, haben ebenfalls einen signifikanten Einfluss auf die Nutzerzahl.
Sieht man jedoch auf den Wert für R-Quadrat, der sich auf 0.3996 beläuft, wird klar, dass das Modell die Situation nicht in ihrer Vollständigkeit erfassen kann. Dieser Wert ist sehr niedrig und sagt aus, dass nur 39.96% der Varianz in den Daten durch dieses Regressionsmodell erklärt werden können. Wirft man einen Blick auf die Visualisierungen im Kapitel Basis Visualisierungen wird klar weshalb. Die lineare Regression kann nur in eine Richtung vorhersagen treffen. Z.B. Wenn sich x um 1 erhöht/verringert, erhöht/verringert sich y um z. Es ist jedoch zu sehen, dass die Zusammenhänge in den Daten nicht linearer Natur sind, denn wenn die Temperatur eine bestimmte Marke überschreitet fahren weniger Menschen mit dem Rad. Das Selbe lässt sich bei der Uhrzeit beobachten. Morgens fahren viele Menschen Rad, über Tag weniger und am Nachmittag wieder mehr, bis es in der Nacht wieder weniger wird. Es ist hier also wenig verwunderlich, dass ein lineares Modell den Einfluss der einzelnen Variablen nicht exakt wiedergibt.
Nun wird mit dem erstellten Modell eine Vorhersage für die Testdaten getroffen und der RMSE gemessen. Diese Regression wird später nicht weiter als Benchmark genutzt, da das LSTM Netz nur für den Bezirk Chelsea erstellt werden kann und eine Vergleichbarkeit so nicht gegeben ist.
#Vorhersage mit Testdaten
lin_reg_q1_nh_predictions <- predict(lin_reg_q1_nh, hourly_test)
linreg_q1_nh_test_mse <- mean((hourly_test$user_count - lin_reg_q1_nh_predictions) ^ 2)
linreg_q1_nh_test_rmse <- sqrt(mean((hourly_test$user_count - lin_reg_q1_nh_predictions) ^ 2))
linreg_q1_nh_test_rmse## [1] 60.80118
Die Vorhersagequalität ist mit einem RMSE von rund 61 Personen pro Stunde und Bezirk nicht sehr hoch. Die Nutzerzahlen bewegen sich überlicherweise in einem niedrigen dreistelligen Bereich. Somit würde eine Planung für die Verteilung der Räder auf die verschiedenen Bezirke, die sich im Schnitt um 61 Personen irrt als nicht gewinnbringend empfunden werden.
Der nachfolgende Code wählt nur die statistisch signifikanten Variablen aus dem vorherigen Modell aus und überführt diese in ein neues Modell. Dies wird im folgenden nicht weiter genutzt.
# Bestimmen welche Variablen signifikant sind
test_sig <- summary(lin_reg_q1_nh)$coefficients[-1,4] < 0.05
# Nur signifikante Behalten
sig_x <- names(test_sig)[test_sig == TRUE]
View(sig_x)
sig_x <- str_replace_all(sig_x, c(" " = "_"))
sig_x <- str_replace_all(sig_x, c("-" = "_"))
sig_x <- str_replace_all(sig_x, c("'" = ""))
# Model mit nur signifikanten Variablen erstellen
sig.model <- as.formula(paste("user_count", paste(sig_x, collapse = " + "), sep = " ~ "))Nachfolgend wir eine lineare Regression auf Stationsebene durchgeführt.
# Regression mit Stationen
model_q1 <- as.formula(user_count ~ .)
lin_reg_q1 <- lm(model_q1, hourly_starts)Da der Datensatz in bearbeiteter Form mehrere Gigabyte groß ist, kann keine Analyse durchgeführt werden. Nachfolgend wird ausschließlich der Datensatz und die Regression für den Bezirk Chelsea betrachtet.
# Regression mit Chelsea
hourly_chelsea_model <- as.formula(user_count ~ start_date + start_hour + quarter + holiday + weekday + temp + wdsp + precip)
hourly_chelsea_lin_reg <- lm(hourly_chelsea_model, hourly_chelsea_train)
summary(hourly_chelsea_lin_reg)##
## Call:
## lm(formula = hourly_chelsea_model, data = hourly_chelsea_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -297.14 -86.47 -18.81 55.31 653.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 412.03227 885.53756 0.465 0.6417
## start_date -0.02342 0.05286 -0.443 0.6577
## start_hour 7.21918 0.20267 35.621 <2e-16 ***
## quarter 8.51357 4.98197 1.709 0.0875 .
## holiday -83.28945 9.32655 -8.930 <2e-16 ***
## weekdayDienstag 7.96353 5.37908 1.480 0.1388
## weekdayMittwoch 11.88078 5.35552 2.218 0.0266 *
## weekdayDonnerstag 10.02969 5.31911 1.886 0.0594 .
## weekdayFreitag -12.09804 5.26397 -2.298 0.0216 *
## weekdaySamstag 6.43192 9.76822 0.658 0.5103
## weekdaySonntag 1.22955 9.78282 0.126 0.9000
## temp 4.71209 0.15037 31.337 <2e-16 ***
## wdsp 0.15433 0.09198 1.678 0.0934 .
## precip -27.93937 3.09001 -9.042 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 129.9 on 8649 degrees of freedom
## Multiple R-squared: 0.2951, Adjusted R-squared: 0.294
## F-statistic: 278.5 on 13 and 8649 DF, p-value: < 2.2e-16
hourly_linreg_train_mse <- mean((hourly_chelsea_lin_reg$residuals)^2)
hourly_linreg_train_rmse <- sqrt(mean((hourly_chelsea_lin_reg$residuals)^2))
hourly_linreg_train_rmse## [1] 129.7962
Die Zusammenfassung des Modells zeigt, dass nur wenige Variablen einen signifikanten Einfluss auf die Nutzerzahl im Bezirk Chelsea zu haben scheinen. Einige Wochentage, das Datum, das Quartal, sowie die Windgeschwindigkeit scheinen hier nur einen zufälligen Einfluss zu haben. Sie werden aus den gleichen Gründen wie bereits bei der Regression mit allen Bezirken dennoch im Modell belassen. Die Güte des Gesamtmodells ist mit einem R-Quadrat von 0.294 als eher schlecht zu bewerten. Nur 29.4% der Varianz in den Daten wird durch das Modell erklärt. Dies ist wie oben bereits erläutert nicht weiter verwunderlich, da viele Zusammenhänge zwischen den Variablen nicht linearer natur zu sein scheinen.
Nun wird mit dem erstellten Modell eine Vorhersage für die Testdaten getroffen und der RMSE gemessen. Diese Regression wird später als Benchmark für das LSTM Netz genutzt.
# Vorhersage Testdaten mit Chelsea
hourly_chelsea_lin_reg_predictions <- predict(hourly_chelsea_lin_reg, hourly_chelsea_test)
hourly_linreg_test_mse <- mean((hourly_chelsea_test$user_count - hourly_chelsea_lin_reg_predictions) ^ 2)
hourly_linreg_test_rmse <- sqrt(mean((hourly_chelsea_test$user_count - hourly_chelsea_lin_reg_predictions) ^ 2))
hourly_linreg_test_rmse## [1] 152.6931
Der RMSE ist mit duchschnittlich rund 153 Personen pro Stunde im Bezirk Chelsea um die sich die Vorhersage irrt sehr hoch. Da auch hier die Nutzeranzahl pro Stunde meist in einem niedrigen dreistelligen Bereich liegt, ist die Vorhersagequalität des linearen Regressionsmodells als schlecht zu bewerten. Besonders deutlich wird dies im visuellen Vergleich. Die Vorhersage nimmt die Schwankungen in den Daten zwar mit, ist aber zu schwach.
hourly_vizdata <- data.frame( Date = hourly_chelsea_test$start_datetime, Values = hourly_chelsea_test$user_count, LinReg = hourly_chelsea_lin_reg_predictions)
ggplot() +
labs(title = "Lineare Regression Vorhersage Stündliche Nutzerzahlen") +
xlab("Date") + ylab("User Count") +
geom_line(data = hourly_vizdata, aes(y = Values, x = Date, color = "black")) +
geom_line(data = hourly_vizdata, aes(y = LinReg, x = Date, color = "orange")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))+
facet_wrap(month(hourly_chelsea_test$start_datetime, label = T, abbr = F), scales = "free_x") +
scale_color_manual(name = "Legende", values = c("black" = "black", "orange" = "orange"), labels = c("Tatsächliche Nutzer", "Lineare Regression Vorhersage"))Tägliche Regression
Hier eine tägliche Reggression für alle Stationen/Neighborhoods.
#Regression mit q2_starts_w
options(max.print=450)
model_q2 <- as.formula(anzahl ~ .-date)
lin_reg_q2 <- lm(model_q2, q2_starts_w)
summary(lin_reg_q2)#Regression mit q2_starts_w2_agg
options(max.print=450)
model_q2_w2_agg <- as.formula(anzahl ~ .-start_date)
lin_reg_q2_w2_agg <- lm(model_q2_w2_agg, q2_starts_w2_agg)
summary(lin_reg_q2_w2_agg)Prediction.
#Vorhersage mit Testdaten
lin_reg_q2_predictions <- predict(lin_reg_q2_w2_agg, test_data_q2)
test_q2_mse <- mean((test_data_q2$anzahl - lin_reg_q2_predictions) ^ 2)
test_q2_rmse <- sqrt(mean((test_data_q2$anzahl - lin_reg_q2_predictions) ^ 2))
test_q2_mse
test_q2_rmseDer RMSE für alle Neighborhoods von der linearen Regression liegt bei 494,8. Das heißt, dass sich die lineare Regression um ca 500 Nutzer pro Neighborhood pro Tag verschätzt. Dies ist in Anbetracht der großen Datenmenge und der vielen Features ein guter Wert. Wenn man von ca. 3000 - 5000 Nutzern pro Station pro Tag ausgeht verschätzt sich die Regression um ca. 10-20%. Dies sind nur geschätzte Werte da die genauen Zahlen sehr variieren.
Tägliche Lineare Regression für ein Neighborhood
Eine lineare Regression für das Neighborhood Chelsea zu den Nutzerzahlen pro Tag.
#Regression mit daily_neighborhood_train
options(max.print=150)
daily_neighborhood_model <- as.formula(anzahl ~ .)
lin_reg_daily_neighborhood <- lm(daily_neighborhood_model, daily_neighborhood_train)
summary(lin_reg_daily_neighborhood)##
## Call:
## lm(formula = daily_neighborhood_model, data = daily_neighborhood_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2008.98 -336.96 82.35 379.06 1810.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.109e+05 6.272e+04 1.767 0.07808 .
## start_date -5.529e+00 3.733e+00 -1.481 0.13951
## avg_tripduration 6.400e-02 1.804e-01 0.355 0.72298
## avg_age -3.689e+02 7.583e+01 -4.865 1.77e-06 ***
## maximum.temperature 5.229e+03 4.580e+03 1.142 0.25434
## minimum.temperature 5.110e+03 4.580e+03 1.116 0.26530
## average.temperature -1.028e+04 9.159e+03 -1.122 0.26263
## precipitation -5.826e+01 4.611e+00 -12.636 < 2e-16 ***
## snow.fall 6.211e+00 6.215e+00 0.999 0.31832
## snow.depth -8.241e+00 1.844e+00 -4.468 1.08e-05 ***
## weekdayDienstag 3.627e+02 1.236e+02 2.935 0.00357 **
## weekdayMittwoch 3.910e+02 1.239e+02 3.155 0.00175 **
## weekdayDonnerstag 3.745e+02 1.236e+02 3.030 0.00264 **
## weekdayFreitag -1.421e+02 1.215e+02 -1.169 0.24315
## weekdaySamstag 7.626e+01 2.264e+02 0.337 0.73648
## weekdaySonntag -2.091e+02 2.298e+02 -0.910 0.36352
## holiday -2.069e+03 2.184e+02 -9.475 < 2e-16 ***
## monthFebruar -2.996e+02 2.137e+02 -1.402 0.16192
## monthMärz 2.952e+01 2.980e+02 0.099 0.92115
## monthApril 1.432e+02 4.032e+02 0.355 0.72277
## monthMai 4.079e+02 5.120e+02 0.797 0.42623
## monthJuni 8.846e+02 6.240e+02 1.418 0.15718
## monthJuli 7.591e+02 7.386e+02 1.028 0.30481
## monthAugust 7.954e+02 8.566e+02 0.929 0.35378
## monthSeptember 1.675e+03 9.599e+02 1.745 0.08187 .
## monthOktober 2.021e+03 1.067e+03 1.895 0.05896 .
## monthNovember 1.675e+03 1.173e+03 1.429 0.15397
## monthDezember 1.564e+03 1.281e+03 1.221 0.22298
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 603.9 on 334 degrees of freedom
## Multiple R-squared: 0.8683, Adjusted R-squared: 0.8576
## F-statistic: 81.53 on 27 and 334 DF, p-value: < 2.2e-16
#Vorhersage mit Testdaten
lin_reg_q2n_predictions <- predict(lin_reg_daily_neighborhood, daily_neighborhood_test)
test_q2_mse <- mean((daily_neighborhood_test$anzahl - lin_reg_q2n_predictions) ^ 2)
test_q2_rmse <- sqrt(mean((daily_neighborhood_test$anzahl - lin_reg_q2n_predictions) ^ 2))
test_q2_mse## [1] 5094109
## [1] 2257.013
Wir bekommen einen RMSE von 2257. D.h., dass die Vorhersage der Regression um ca 2257 Nutzer daneben liegt pro Tag. Dies ist ein sehr schlechtes Ergebnis. Man könnte jetzt Features entfernen die nicht statistisch signifikant sind und die Regression erneut durchlaufen lassen, um zu gucken, ob dies zu einer Verbesserung führt. Ich überspringe diesen Punkt und gehe weiter zur Visualisierung und zum LSTM.
daily_vizdata <- data.frame( Date = daily_neighborhood_test$start_date, Values = daily_neighborhood_test$anzahl, LinReg = lin_reg_q2n_predictions)
ggplot() +
labs(title = "Lineare Regression Vorhersage Täglicher Nutzerzahlen") +
xlab("Date") + ylab("Anzahl") +
geom_line(data = daily_vizdata, aes(y = Values, x = Date, color = "black")) +
geom_line(data = daily_vizdata, aes(y = LinReg, x = Date, color = "orange")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))+
facet_wrap(month(daily_neighborhood_test$start_date, label = T, abbr = F), scales = "free_x") +
scale_color_manual(name = "Legende", values = c("black" = "black", "orange" = "orange"), labels = c("Tatsächliche Nutzer", "Lineare Regression Vorhersage"))Man kann in der Vorhersage sehen das die Schwankungen von der Regression erkannt werden, aber die Werte leider zu stark abweichen. Wie oben schon erwähnt, müsste man diese Regression erst optimieren, damit sie aussagekräftig ist.
Künstliches Neuronales Netz
Zur Vorhersage der Nutzerzahlen pro Tag und pro Stunde wird eine spezielle Form der künstlichen neuronalen Netze genutzt, das Long Short Term Memory Netz. Diese Sonderform ist besonders für das Arbeiten mit Zeitreihen gerüstet, da die einzelnen Neuronen in den Schichten des Netzes vorherige Zustände speichern und so mit in die Berechnung der einzelnen Gewichte einfließen lassen können.
Sowohl die täglichen Daten, als auch die stündlichen Daten stellen eine multivariate Zeitreihe dar. Dies bedeutet, dass eine Größe y, hier die Nutzeranzahl, im Zeitablauf, hier Tage bzw. Stunden, nicht nur von der Zeit an sich abhängt, sondern zusätzlich auch von weiteren Features x, hier z.B. das Quartal oder die Temperatur.
Um gute Vorhersageergebnisse zu erzielen, sind hier besondere Schritte notwendig, bevor die Daten in das LSTM Netz eingespeist werden.
Zunächst wird die Größe y, hier die Nutzeranzahl, stationarisiert. Das bedeutet, dass anstelle der realen Werte die Differenzen von einem Zeitpunkt zum nächsten betrachtet werden, also der Nutzerzuwachs oder -verlust. Dies wird gemacht und schein-korrelativen und schein-regressiven Effekten vorzubeugen und Trends und Saisonalitäten aus der Zeitreihe zu eliminieren.
Außerdem werden Lags für die weiteren Features x erstellt. Dies bedeutet, dass alle Werte einer Periode t eine Periode nach hinten in die Periode t+1 verschoben werden. Alle features aus der Periode t, sowie die gelaggten features werden in den Berechnungen berücksichtigt. Dies erhöht die Anzahl der Features, die die Größe y beeinflussen. Wie viele lags pro Feature optimal sind kann mittels der Methode ARIMAX ermittelt werden. Es kann außerdem überlegt werden, ob für die differenzierten y-Werte ebenfalls Lags gebildet werden, die dann als Features in die Periode t aufgenommen werden. Dies kann die Prognosegüte erhöhen, hat allerdings zur Folge, dass mehr Fälle ausgeschlossen werden, da sich durch die Verschiebung fehlende Werte in den ersten Zeilen ergeben. Aus der Anzahl der Lags ergibt sich die Anzahl an Timesteps, die für das LSTM Netz eine wichtige Komponente ist. Im multivariaten Fall ergibt sich die Anzahl der Timesteps aus der Anzahl der Lags+1. Jedes Lag mit allen Variablen darin ergibt einen timestep. Die Periode t wird ebenfalls dazu gezählt.
Die umstrukturierten Daten werden dann in die Form eines dreidimensionalen Arrays überführt, da LSTMs diese als Input erwarten. Das dann trainierte Modell sagt danach nicht die tatsächlichen y-Werte der Testdaten voraus, sondern die Differenzwerte von einer Periode zur nächsten. Dies muss am Ende umgerechnet werden, um ein Ergebnis zu erhalten, das leichter zu interpretieren und nützlich ist.
Stündliche Daten Aufbereiten
Um ein künstliches neuornales Netz zu erhalten, das fehlerfrei läuft, müssen der Trainings- und Testdatensatz der stündlichen Fahrten noch weiter aufbereitet werden. Für die Spalten neighborhood, weekday, month, week, quarter und start_hour werden Dummyvariablen erstellt. Dies ist notwendig, da es sich um kategoriale Variablen handelt. Da eine reine Quantifizerung der Variablen durch hohe Zahlenwerte für einige Kategorein eine Verzerrung der Analysen zur Folge hätte, werden binäre Variblen für jede Kategorie erstellt, die die Werte 0 und 1 annehmen. Um bei der späteren Anwendung der künstlichen neuronalen Netze keine Fehler zu produzieren, werden die Dummys effektkodiert. Statt des Wertes 0 wird der Wert -1 genutzt. Dieser Vorgang vergrößert die Datensätze enorm, da für jede Kategorie eine Dummyvariable erstellt werden muss, z.B. 23 Dummyvariablen für das Merkmal Stunde.
Dies wird hier für die gesamten Trainings- und Testdatensätze, sowie für die Trainings- und Testdatensätze, die nur den Bezirk Chelsea betreffen, durchgeführt. Nur die Datensätze für den Bezirk Chelsea werden weitergehen betrachtet und für die Analyse genutzt.
# Dummyvariablen erstellen
hourly_starts_nh <- dummy_cols(hourly_starts_nh, c("start_hour", "neighborhood", "weekday", "month", "week"), remove_first_dummy = T, remove_selected_columns = T)
hourly_test <- dummy_cols(hourly_test, c("start_hour", "neighborhood", "weekday", "month", "week"), remove_first_dummy = T, remove_selected_columns = T)
# Namen der Dummys bearbeiten
names(hourly_starts_nh) <- str_replace_all(names(hourly_starts_nh), c(" " = "_"))
names(hourly_test) <- str_replace_all(names(hourly_test), c(" " = "_"))
# Dummys effektkodieren
hourly_starts_nh[,c(3, 7:146)] <- effectcoding(hourly_starts_nh[,c(3, 7:146)])
hourly_test[,c(3, 7:105)] <- effectcoding(hourly_test[,c(3, 7:105)])# Chelsea Trainingsdaten dummyfizieren und effektcodieren
hourly_chelsea_train <- dummy_cols(hourly_chelsea_train, c("start_hour", "weekday", "quarter"), remove_first_dummy = T, remove_selected_columns = T)
names(hourly_chelsea_train) <- str_replace_all(names(hourly_chelsea_train), c(" " = "_"))
hourly_chelsea_train[,c(5, 9:40)] <- effectcoding(hourly_chelsea_train[,c(5, 9:40)])
# Chelsea Testdaten dummyfizieren und effektcodieren
hourly_chelsea_test <- dummy_cols(hourly_chelsea_test, c("start_hour", "weekday", "quarter"), remove_first_dummy = T, remove_selected_columns = T)
names(hourly_chelsea_test) <- str_replace_all(names(hourly_chelsea_test), c(" " = "_"))
hourly_chelsea_test[,c(5, 9:40)] <- effectcoding(hourly_chelsea_test[,c(5, 9:40)])ARIMAX
Um für die Trainingsdaten die optimale Anzahl an Lags zu bestimmen, werden die vorbereiteten Trainingsdaten der auto.arima Funktion übergeben, die mit Hilfe der ARIMAX Methode die optimale Anzahl an Lags ermittelt. Daraus wird dann die Anzahl der Timesteps generiert.
# ARIMAX für Bestimmung der lags mit Trainingsdaten
hourly.arimax.Y <- as.numeric(hourly_chelsea_train$user_count)
hourly.arimax.X <- as.matrix(hourly_chelsea_train[,c(5:40)])
hourly_arimax <- auto.arima(hourly.arimax.Y, xreg = hourly.arimax.X)
hourly_lags <- hourly_arimax$arma[1]
hourly.ann.lags <- as.lags(hourly_lags, type = "multivariate")
hourly.ann.timesteps <- as.timesteps(hourly.ann.lags, type = "multivariate")Die optimale Laganzahl hier ist 5. Daraus ergibt sich eine Anzahl von Timesteps von 6. Diese gehen als Hyperparameter in das LSTM-Modell mit ein.
Resampling Und Normalisieren Der Trainings- und Testdaten
Nun wird die Anzahl der Radfahrer stationarisiert, sowie Lags für die Differnzreihe gebildet, die als weitere Features mit in das LSTM Modell eingehen sollen. Diese werden zu den Trainingsdaten hinzugefügt. Danach wird die y-Variable user_count_diff, sowie die Features getrennt voneinander geresampled, sodass die Lags in den Datensatz mitaufgenommen werden. Danach können sie wieder in einem Data Frame zusammengefasst werden. Danach wird das Datum und die Uhrzeit wieder hinzugefügt, um einen besseren Überblick über den Datensatz zu erhalten. Dies wird analog auch für die Testdaten durchgeführt.
# Trainingsdaten mit gelaggtem y als x vorbereiten
hourly_tsdiff <- as.data.frame(build_stationary(hourly_chelsea_train$user_count, k = hourly.ann.lags))
if (hourly.ann.lags > 0) {
hourly_tsdiff <- hourly_tsdiff[c(-1:-hourly.ann.lags),]
colnames(hourly_tsdiff) <- c("user_count_diff", "user_count_diff_lag1","user_count_diff_lag2",
"user_count_diff_lag3", "user_count_diff_lag4","user_count_diff_lag5")
}
hourly.df.ann <- cbind.data.frame(hourly_chelsea_train[c((hourly.ann.lags+2):nrow(hourly_chelsea_train)),], hourly_tsdiff)
hourly.df.ann <- hourly.df.ann[,c(1,3:4,41:46,5:40)]
hourly.Y.train.resample <- resample.y(hourly.df.ann$user_count_diff, hourly.ann.timesteps)
hourly.X.train.resample <- resample.X(hourly.df.ann[,5:45], hourly.ann.timesteps)
hourly.df.resample <- cbind.data.frame(hourly.Y.train.resample, hourly.X.train.resample)
colnames(hourly.df.resample)[1] <- c("user_count_diff")
hourly.df.resample$start_datetime <- hourly.df.ann$start_datetime[(hourly.ann.timesteps):nrow(hourly.df.ann)]
hourly.df.resample <- hourly.df.resample[,c(248, 1:247)]
# Testdaten mit gelaggtem y als x vorbereiten
hourly_tsdiff_test <- as.data.frame(build_stationary(hourly_chelsea_test$user_count, k = hourly.ann.lags))
if (hourly.ann.lags > 0) {
hourly_tsdiff_test <- hourly_tsdiff_test[c(-1:-hourly.ann.lags),]
colnames(hourly_tsdiff_test) <- c("user_count_diff", "user_count_diff_lag1","user_count_diff_lag2",
"user_count_diff_lag3", "user_count_diff_lag4","user_count_diff_lag5")
}
hourly.df.ann.test <- cbind.data.frame(hourly_chelsea_test[c((hourly.ann.lags+2):nrow(hourly_chelsea_test)),], hourly_tsdiff_test)
hourly.df.ann.test <- hourly.df.ann.test[,c(1,3:4,41:46,5:40)]
hourly.Y.test.resample <- resample.y(hourly.df.ann.test$user_count_diff, hourly.ann.timesteps)
hourly.X.test.resample <- resample.X(hourly.df.ann.test[,5:45], hourly.ann.timesteps)
hourly.df.resample.test <- cbind.data.frame(hourly.Y.test.resample, hourly.X.test.resample)
colnames(hourly.df.resample.test)[1] <- c("user_count_diff")
hourly.df.resample.test$start_datetime <- hourly.df.ann.test$start_datetime[(hourly.ann.timesteps):nrow(hourly.df.ann.test)]
hourly.df.resample.test <- hourly.df.resample.test[,c(248, 1:247)]
hourly.df.train <- hourly.df.resample
hourly.df.test <- hourly.df.resample.testUm nach der Prognose der Differenzen der Nutzerzahlen die absoluten Zahlen wiederherstellen zu können, werden die originären y-Werte der Testdaten gespeichert.
# Originäre y-Werte der Testdaten speichern für späteres re-resampling
hourly_y_test <- hourly_chelsea_test$user_count[12:NROW(hourly_chelsea_test)]
hourly_y_prior_first_test <- hourly_chelsea_test$user_count[11]Im nächsten Schritt werden alle Variablen beider Datensätze, die keine Dummyvariablen sind, normalisiert, damit keine unproportionalen Einflüsse einiger Variablen entstehen.
# Trainings- und Testdaten normalisieren
hourly_selected_columns <- which(colnames(hourly.df.train)%in%c("user_count_diff",
"temp_1","wdsp_1","precip_1","user_count_diff_lag1_1","user_count_diff_lag2_1","user_count_diff_lag3_1","user_count_diff_lag4_1","user_count_diff_lag5_1", "temp_2","wdsp_2","precip_2","user_count_diff_lag1_2","user_count_diff_lag2_2","user_count_diff_lag3_2","user_count_diff_lag4_2","user_count_diff_lag5_2", "temp_3","wdsp_3","precip_3","user_count_diff_lag1_3","user_count_diff_lag2_3","user_count_diff_lag3_3","user_count_diff_lag4_3","user_count_diff_lag5_3", "temp_4","wdsp_4","precip_4","user_count_diff_lag1_4","user_count_diff_lag2_4","user_count_diff_lag3_4","user_count_diff_lag4_4","user_count_diff_lag5_4", "temp_5","wdsp_5","precip_5","user_count_diff_lag1_5","user_count_diff_lag2_5","user_count_diff_lag3_5","user_count_diff_lag4_5","user_count_diff_lag5_5", "temp_6","wdsp_6","precip_6","user_count_diff_lag1_6","user_count_diff_lag2_6","user_count_diff_lag3_6","user_count_diff_lag4_6","user_count_diff_lag5_6"))
hourly_nd <- normalize_data(hourly.df.train[,hourly_selected_columns], hourly.df.test[,hourly_selected_columns])
hourly.df.train.norm <- hourly.df.train
hourly.df.train.norm[,hourly_selected_columns] <- hourly_nd$train
hourly.df.train.norm <- hourly.df.train.norm[,-1]
hourly.df.test.norm <- hourly.df.test
hourly.df.test.norm[,hourly_selected_columns] <- hourly_nd$test
hourly.df.test.norm <- hourly.df.test.norm[,-1]Sind die Daten normalisiert ist die Vorbereitung abgeschlossen und es können die Matrizen als Input für das neuronale Netz erstellt werden. X- und Y-Variablen werden hierbei getrennt.
Tägliche Daten Aufbereiten Für Alle Neighborhoods
Vorbereiten des Datensatzes für ein LSTM.
# Data Table variante
#library(data.table)
#q2_starts_w_nn <- as.data.table(q2_starts_w)
# Dataframe für LSTM erstellen
q2_starts_w_nn <- q2_starts_w
#Dummyvariablen erstellen
q2_starts_w_nn <- dummy_cols(q2_starts_w, c("weekday", "start_station_name"), remove_first_dummy = T, remove_selected_columns = T)
#Namen der Dummys bearbeiten
names(q2_starts_w_nn) <- str_replace_all(names(q2_starts_w_nn), c(" " = "_"))
#Dummys effektkodieren
#q2_starts_w_nn[,c(11:670)] <- effectcoding(q2_starts_w_nn[,c(11:670)])Funktioniert nicht! Da zuviele Features das RStudio zum Abstürzen bringt.
Neue Variante mit Neighborhoods.
Trainingdaten Für LSTM Vorbereiten
Dummyvariablen erstellen und bearbeiten, um sie dann zu effektcodieren.
# Data Table variante
#library(data.table)
#q2_starts_w2_nn <- as.data.table(q2_starts_w2_agg)
# Dataframe für LSTM erstellen
q2_starts_w2_nn <- q2_starts_w2_agg
#Dummyvariablen erstellen
q2_starts_w2_nn <- dummy_cols(q2_starts_w2_nn, c("weekday", "neighborhood", "month"), remove_first_dummy = T, remove_selected_columns = T)
#Namen der Dummys bearbeiten
names(q2_starts_w2_nn) <- str_replace_all(names(q2_starts_w2_nn), c(" " = "_"))
#Dummys effektkodieren
q2_starts_w2_nn[,c(11:76)] <- effectcoding(q2_starts_w2_nn[,c(11:76)])Testdaten Für LSTM Vorbereiten
Gleiche Schritte wie bei den Trainingsdaten durchführen.
# Data Table variante
#test_data_q2_nn <- as.data.table(test_data_q2)
# Dataframe für LSTM erstellen
test_data_q2_nn <- test_data_q2
#Dummyvariablen erstellen
test_data_q2_nn <- dummy_cols(test_data_q2_nn, c("weekday", "neighborhood", "month"), remove_first_dummy = T, remove_selected_columns = T)
#Namen der Dummys bearbeiten
names(test_data_q2_nn) <- str_replace_all(names(test_data_q2_nn), c(" " = "_"))
#Dummys effektkodieren
test_data_q2_nn[,c(11:68)] <- effectcoding(test_data_q2_nn[,c(11:68)])Vergleich der Datensätze.
## character(0)
test_data_q2_nn$neighborhood_Astoria <- NULL
test_data_q2_nn$neighborhood_Crown_Heights <- NULL
test_data_q2_nn$neighborhood_Ditmars_Steinway <- NULL
test_data_q2_nn$neighborhood_Harlem <- NULL
test_data_q2_nn$neighborhood_Morningside_Heights <- NULL
test_data_q2_nn$neighborhood_Prospect_Heights <- NULL
test_data_q2_nn$'neighborhood_Prospect-Lefferts_Gardens' <- NULL
rm(df1,df2)ARIMAX Für Alle Neighborhoods Tägliche Nutzer
Aufwand zu groß alle Neighborhoods in einem Netz vorauszusagen.
# Build ARIMAX model for training Data
require(forecast)
arx.Y <- q2_starts_w2_nn$anzahl
arx.X <- as.matrix(q2_starts_w2_nn[,c(3:76)])
arx <- auto.arima(arx.Y, xreg = arx.X)
arx
non_seasonal_ar_order <- arx$arma[1] # p
non_seasonal_ma_order <- arx$arma[2] # q
seasonal_ar_order <- arx$arma[3]
seasonal_ma_order <- arx$arma[4]
period_of_data <- arx$arma[5]
non_seasonal_diff_order <- arx$arma[6] #d
seasonal_diff_order <- arx$arma[7]LSTM Für Alle Neighborhoods Tägliche Nutzer
######
# LSTM
# Hyperparameters
ann.epochs <- 100
ann.batchsize <- 32
ann.lags <- as.lags(non_seasonal_ar_order, type = "multivariate")
ann.timesteps <- as.timesteps(ann.lags, type = "multivariate")
# Transform y (= price) into stationary and create lagged dataset
tsdiff <- as.data.frame(build_stationary(q2_starts_w2_nn$anzahl, k = ann.lags))
'if (ann.lags > 0) {
tsdiff <- tsdiff[c(-1:-ann.lags),]
colnames(tsdiff) <- c("price_diff", "pdiff_lag1","pdiff_lag2","pdiff_lag3","pdiff_lag4")
}'
df.ann <- cbind.data.frame(q2_starts_w2_nn[c((ann.lags+2):nrow(q2_starts_w2_nn)),], tsdiff)
df.ann <- df.ann[,c(1:4,45:49,5:44)]Uns ist beim kreieren des Datensatzes für das LSTM aufgefallen, dass wir durch die vielen Neighborhoods ein Problem kriegen. Deswegen sind wir zu dem Punkt gekommen, erstmal ein LSTM für ein Neighborhood exemplarisch zu bauen.
Long Short Term Memory Netz für stündliche Daten Erstellen
Bevor das Modell endgültig erstellt und trainiert werden kann, werden die letzten Hyperparameter festgelegt. Nach sehr vielen Versuchen haben eine Batchsize von 1 und eine Epochenanzahl von 120 das beste Ergebnis erzielt.
Die Funktion fit_lstm erstellt nun das LSTM Modell und trainiert es anschließend mit den gesetzten Parametern. Als Inputdaten werden die X und Y Matrizen eingesetzt. Diese werden automatisch von der Funktion in dreidimensionale Arrays umgewandelt, die als Input erwartet werden.
Die Anzahl der Timesteps ist in diesem Fall 6 und wurde in der unten angegebenen Variable gespeichert. Auch die Anzahl der Epochen und die Batchsize wird angegeben. Die Batchsize wird unterteil in den Wert NA und die selbst festgelegte Größe. Dies ist der Fall, da das Netz hier auch in der Inputschicht eine Batchsize erwartet. Diese muss im Falle einer multivariaten Zeitreihe nicht zwingend angegeben werden und wird daher mit NA angegeben.
Die nächsten Argumente können genutzt werden, um mit Hilfe der k-fold cross validation die optimale Anzahl an Epochen zu bestimmen. In diesem Fall hat die Nutzung dieser Methodik zu schlechteren Ergebnissen geführt als die Epochenanzahl von 120, weshalb diese beibehalten wurde.
Im nächsten Argument hidden werden die Anzahl, Größe und Aktivierungsfunktion der Hiddenschichten spezifiziert. In diesem Fall hat nach vielen erfolglosen Versuchen das Netz mit zwei Hiddenschichten mit jeweils 50 und 25 Neuronen, sowie der Tangens-Hyperbolicus Aktivierungsfunktion die besten Ergebnisse erzielt.
Um Overfitting zu vermeiden und diesem entgegenzuwirken, werden zwei sog. Dropout-Layer eingefügt. Diese befinden sich zwischen den einzelnen Schichten und filtern nach dem Zufallsprinzip einen bestimmten Prozentsatz Inputneuronen heraus, die den Wert 0 erhalten. Hier wurden mit einer Dropout-Rate von jeweils 50% und 10% die besten Ergebnisse erzielt.
Die Aktivierungsfunktion für den Output-Layer ist linear gewählt, da es sich um ein Regressionsproblem handelt, nicht um ein Klassifizierungsproblem. Es wird kein stateful processing genutzt, da hierbei in den einzelnen Neuronen ausschließlich der letzte Wert genutzt wird. Dies empfiehlt sich für univariate Zeitreihen, nicht aber für multivariate. Das Argument return_sequences legt hier fest, dass Vorhersagewerte für alle Timesteps erstellt werden.
Als Maßeinheit für Fehler wird der Mean Squared Error (MSE) gewählt. Zusätzlich wird durch das Argument metrics eine zweite Maßeinheit gewählt. Hier ist dies der Mean Absolute Error (MAE).
Der Optimierungsalgorithmus legt fest wie genau das LSTM Netz lernt. Hier hat der Follow-The-Regularized-Leader, kurz FTRL, Algorithmus die mit Abstand besten Ergebnisse geliefert. Dieser wird in der Praxis häufig für large-scale Probleme verwendet, wie z.B. das Vorhersagen von Klickzahlen auf Werbebanner.
# Stündliches LSTM trainieren
hourly_lstm <- fit_lstm(X = hourly.X.train, y = hourly.Y.train,
timesteps = hourly.ann.timesteps,
validation_split = 0.05,
epochs = hourly.ann.epochs,
batchsize = c(NA,hourly.ann.batchsize),
#k.fold = 3, k.optimizer = "min",
hidden = data.frame(c(50,25),c("tanh")),
dropout = c(0.5,0.1),
output_activation = "linear",
stateful = F, return_sequences = T,
loss = "mean_squared_error",
optimizer = "Ftrl",
metrics = c('mean_absolute_error'))Nun, da das Netz trainiert ist, können die Nutzerzahlen der Tesdatenmenge vorhergesagt werden. Hier werden zunächst die Differenzen vorhergesagt und im Anschluss denormalisiert. Danach wird mit Hilfe der gespeicherten originalen Nutzerzahlen die Differenzierung invertiert, sodass statt der Differenzen die “tatsächlichen” Werte betrachtet werden. Danach wird wie bei der linearen Regression der RMSE berechnet.
# Nutzerdifferenzen vorhersagen
hourly.lstm.X.test <- as.LSTM.X(hourly.X.test, hourly.ann.timesteps)
hourly.predictiondiff <- hourly_lstm$model %>% predict(hourly.lstm.X.test, batch_size = hourly.ann.batchsize)
hourly.predictiondiff <- denormalize(hourly.predictiondiff[,1], hourly_nd$min[1], hourly_nd$max[1])
# Differenzierung rückgängig machen
hourly_testseries <- c(hourly_y_prior_first_test, hourly_y_test)
hourly_testseries <- hourly_testseries[-NROW(hourly_testseries)]
hourly_predictions <- invert_differencing(hourly.predictiondiff, hourly_testseries)
# RMSE
hourly_lstm_rmse <- rmse(hourly_testseries, hourly_predictions)
hourly_lstm_rmse## [1] 4.25011
Der RMSE verändert sich geringfügig mit jedem Mal, das das Netz lernt, bewegt sich aber in einem Bereich von rund 2 bis 6 Nutzern, die im Schnitt pro Stunde falsch vorhergesagt werden. Dies ist ein sehr gutes Ergebnis und wird im nächsten Schritt visualisiert.
# Visualisierung der Ergebnisse
hourly_graphdata <- data.frame(Date = hourly.df.test$start_datetime, Value = hourly_testseries, Predictions = hourly_predictions, LinReg = hourly_chelsea_lin_reg_predictions[12:NROW(hourly_chelsea_lin_reg_predictions)])
ggplot() +
labs(title = "LSTM Vorhersage Stündliche Nutzerzahlen") +
xlab("Date") + ylab("User Count") +
geom_line(data = hourly_graphdata, aes(y = Value, x = Date, color = "black")) +
geom_line(data = hourly_graphdata, aes(y = Predictions, x = Date, color = "cadetblue")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))+
facet_wrap(month(hourly.df.test$start_datetime, label = T, abbr = F), scales = "free_x") +
scale_color_manual(name = "Legende", values = c("black" = "black", "cadetblue" = "cadetblue"), labels = c("Tatsächliche Nutzer", "LSTM Vorhersage"))Die Grafik zeigt die stundengenaue Vorhersage der Nutzerzahlen je Monat im Testdatensatz. Die schwarze Linie zeigt die tatsächlichen Nutzerzahlen, die blaue Linie die Vorhergesagten. Es ist zu sehen, dass beide Linien fast deckungsgleich verlaufen. Teilweise weichen die prognostizierten Werte geringfügig von den Spitzen ab. Das Ergebnis ist jedoch als sehr gut zu beurteilen.
Long Short Term Memory Netz Für Tägliche Daten Für Ein Neighborhood
Traininsdaten Für Das LSTM Vorbereiten
Es wird ein Dataframe für das LSTM erstellt. Dieser wir dummyfiziert und effektcodiert.
# Dataframe für LSTM erstellen
daily_neighborhood_train_nn <- daily_neighborhood_train
#Dummyvariablen erstellen
daily_neighborhood_train_nn <- dummy_cols(daily_neighborhood_train_nn, c("weekday", "month"), remove_first_dummy = T, remove_selected_columns = T)
#Namen der Dummys bearbeiten
names(daily_neighborhood_train_nn) <- str_replace_all(names(daily_neighborhood_train_nn), c(" " = "_"))
#Dummys effektkodieren
daily_neighborhood_train_nn[,c(11:28)] <- effectcoding(daily_neighborhood_train_nn[,c(11:28)])Testdaten Für Das LSTM Vorbereiten
Das Gleiche einmal für die Testdaten. Da die Datensätze von den Features gleich sein müssen, um später im LSTM keine Fehler zu verursachen, füge ich die fehlenden Variablen ein damit sie später in der Dummyfizierung als Spalten auftreten.
# Dataframe für LSTM erstellen
daily_neighborhood_test_nn <- daily_neighborhood_test
# Wegen Gleichheit tricksen
df_addrows <- data.frame(start_date= as_date(c("2017-01-01", "2017-03-01", "2017-04-01", "2017-06-01", "2017-07-01", "2017-09-01", "2017-10-01", "2017-12-01")), anzahl = c(NA,NA,NA,NA,NA,NA,NA,NA), avg_tripduration = c(NA,NA,NA,NA,NA,NA,NA,NA), avg_age = c(NA,NA,NA,NA,NA,NA,NA,NA), maximum.temperature = c(NA,NA,NA,NA,NA,NA,NA,NA), minimum.temperature = c(NA,NA,NA,NA,NA,NA,NA,NA), average.temperature = c(NA,NA,NA,NA,NA,NA,NA,NA), precipitation = c(NA,NA,NA,NA,NA,NA,NA,NA), snow.fall = c(NA,NA,NA,NA,NA,NA,NA,NA), snow.depth = c(NA,NA,NA,NA,NA,NA,NA,NA), weekday = c("Montag","Donnerstag","Donnerstag","Donnerstag","Donnerstag","Donnerstag","Donnerstag","Donnerstag"), holiday = c(NA,NA,NA,NA,NA,NA,NA,NA), month = c("Januar", "März", "April", "Juni", "Juli", "September", "Oktober", "Dezember"))
daily_neighborhood_test_nn <- rbind(daily_neighborhood_test_nn, df_addrows)
daily_neighborhood_test_nn <- daily_neighborhood_test_nn %>% arrange(start_date)
#Dummyvariablen erstellen
daily_neighborhood_test_nn <- dummy_cols(daily_neighborhood_test_nn, c("weekday", "month"), remove_first_dummy = F, remove_selected_columns = T)
daily_neighborhood_test_nn <- na.omit(daily_neighborhood_test_nn)
daily_neighborhood_test_nn$month_Januar <- NULL
daily_neighborhood_test_nn$weekday_Montag <- NULL
daily_neighborhood_test_nn <- daily_neighborhood_test_nn[,c(1:12,15,13,14,16:17,21,25,18,24,23,22,19,28,27,26,20 ) ]
#Namen der Dummys bearbeiten
names(daily_neighborhood_test_nn) <- str_replace_all(names(daily_neighborhood_test_nn), c(" " = "_"))
#Dummys effektkodieren
daily_neighborhood_test_nn[,c(11:28)] <- effectcoding(daily_neighborhood_test_nn[,c(11:28)])ARIMAX Für Ein Neighborhood
Mit Arimax findet man die optimalen Lags für ein LSTM Netz. Es ergibt sich ein optimaler Lag von 1.
# Baue ARIMAX Model für Trainingsdaten
require(forecast)
n.arx.Y <- daily_neighborhood_train_nn$anzahl
n.arx.X <- as.matrix(daily_neighborhood_train_nn[,c(3:28)])
n.arx <- auto.arima(n.arx.Y, xreg = n.arx.X)
n.arx## Series: n.arx.Y
## Regression with ARIMA(1,0,0) errors
##
## Coefficients:
## ar1 intercept avg_tripduration avg_age maximum.temperature
## 0.4865 9412.934 0.0741 -149.4483 4359.285
## s.e. 0.0534 2802.354 0.1389 75.5417 3468.148
## minimum.temperature average.temperature precipitation snow.fall
## 4256.656 -8556.985 -54.9977 3.7631
## s.e. 3469.204 6937.222 3.6057 4.7597
## snow.depth holiday weekday_Dienstag weekday_Mittwoch
## -10.0403 -752.2029 191.7117 204.2887
## s.e. 2.1432 94.5897 44.2931 53.0712
## weekday_Donnerstag weekday_Freitag weekday_Samstag weekday_Sonntag
## 189.9768 -57.1360 -149.7586 -250.9964
## s.e. 56.0400 55.2564 93.1457 87.5545
## month_Februar month_März month_April month_Mai month_Juni month_Juli
## -237.0331 -82.3370 -31.3381 5.4641 222.2862 126.3858
## s.e. 131.6358 137.8449 148.5141 156.4829 171.4947 184.7374
## month_August month_September month_Oktober month_November
## 130.3107 347.5635 535.7486 170.1526
## s.e. 194.0376 179.8655 169.8022 155.1303
## month_Dezember
## -30.7913
## s.e. 137.8689
##
## sigma^2 estimated as 297948: log likelihood=-2780.67
## AIC=5619.33 AICc=5624.57 BIC=5732.19
LSTM Für Ein Neighborhood
Nun kommen wir zum Herzstück unserer Forschungsfrage, dem LSTM. Als erstes setze ich hier die Hyperparameter fest. Ich habe mich nach etlichen Probeläufen für eine Epochenanzahl von 200 und eine Batchsize von 2 entschieden. Der Lag und der Timestep sind vorgegeben durch ARIMAX. Ich nehme die Nutzeranzahl und baue daraus eine stationäre Reihe und füge dieser einen Lag hinzu.
######
# LSTM
# Hyperparameters
ann.epochs <- 200
ann.batchsize <- 2
ann.lags <- as.lags(non_seasonal_ar_order, type = "multivariate")
ann.timesteps <- as.timesteps(ann.lags, type = "multivariate")
# Transform y (= anzahl) into stationary and create lagged dataset
tsdiff <- as.data.frame(build_stationary(daily_neighborhood_train_nn$anzahl, k = ann.lags))
if (ann.lags > 0) {
tsdiff <- tsdiff[c(-1:-ann.lags),]
colnames(tsdiff) <- c("anzahl_diff", "adiff_lag1")
}
df.ann <- cbind.data.frame(daily_neighborhood_train_nn[c((ann.lags+2):nrow(daily_neighborhood_train_nn)),], tsdiff)
df.ann <- df.ann[,c(1:2,29:30,3:28)]Resamplen der stationären Daten und hinzufügen von Datuminformation für die spätere Auswertung.
# Neusortieren von stationären Daten weil die Variablen exogen sein müssen.
Y.resample <- resample.y(df.ann$anzahl_diff, ann.timesteps)
X.resample <- resample.X(df.ann[,4:30], ann.timesteps)
df.resample <- cbind.data.frame(Y.resample, X.resample)
colnames(df.resample)[1] <- c("anzahl_diff")
# Erweitere Daten um Datumsinformation für späütere Auswertung
df.resample$datehour <- df.ann$start_date[(ann.timesteps):nrow(df.ann)]
df.resample <- df.resample[,c(56, 1:55)]
df.train <- df.resampleEinmal die stationären Daten einfügen und das Resamplen für die Testdaten durchführen.
tsdiff_t <- as.data.frame(build_stationary(daily_neighborhood_test_nn$anzahl, k = ann.lags))
if (ann.lags > 0) {
tsdiff_t <- tsdiff_t[c(-1:-ann.lags),]
colnames(tsdiff_t) <- c("anzahl_diff", "adiff_lag1")
}
df.ann_t <- cbind.data.frame(daily_neighborhood_test_nn[c((ann.lags+2):nrow(daily_neighborhood_test_nn)),], tsdiff_t)
df.ann_t <- df.ann_t[,c(1:2,29:30,3:28)]
# Neusortieren von stationären Daten weil die Variablen exogen sein müssen.
Y.resample_t <- resample.y(df.ann_t$anzahl_diff, ann.timesteps)
X.resample_t <- resample.X(df.ann_t[,4:30], ann.timesteps)
df.resample_t <- cbind.data.frame(Y.resample_t, X.resample_t)
colnames(df.resample_t)[1] <- c("anzahl_diff")
# Erweitere Daten um Datumsinformation für späütere Auswertung
df.resample_t$datehour <- df.ann_t$start_date[(ann.timesteps):nrow(df.ann_t)]
df.resample_t <- df.resample_t[,c(56, 1:55)]
df.test <- df.resample_tNachdem jetzt die Trainingsdaten und die Testdaten fertig sind, kann ich die Originaldaten speichern um später diese zu nutzen, um zu redifferenzieren. Danach normalisiere ich die Daten in den ausgewählten Spalten. Danach extrahiere ich die X- und die Y-Daten jeweils für die Trainings- und Testdaten.
test.start <- as.POSIXct("2017-02-04")
# Speichere original y Werte und den vorherrigen Wert des ersten Testwertes
y_test <- daily_neighborhood_test_nn$anzahl[(daily_neighborhood_test_nn$start_date >= test.start)] # original Bike Werte
y_prior_first_test <- daily_neighborhood_test_nn$anzahl[which(daily_neighborhood_test_nn$start_date == test.start)-1] # für das redifferencing
# Normalisiere trainings und test Datensets bis auf das Datum und die Dummy-Features
selected_columns <- which(colnames(df.train)%in%c("anzahl_diff", "adiff_lag1_1", "avg_tripduration_1", "avg_age_1",
"maximum.temperature_1", "minimum.temperature_1", "average.temperature_1",
"precipitation_1", "snow.fall_1", "snow.depth_1",
"adiff_lag1_2", "avg_tripduration_2", "avg_age_2", "maximum.temperature_2",
"minimum.temperature_2", "average.temperature_2", "precipitation_2", "snow.fall_2",
"snow.depth_2"
))
nd <- normalize_data(df.train[,selected_columns], df.test[,selected_columns]) # c(2:11,22:30)
df.train.norm <- df.train
df.train.norm[, selected_columns] <- nd$train
df.train.norm <- df.train.norm[,-1]
df.test.norm <- df.test
df.test.norm[, selected_columns] <- nd$test
df.test.norm <- df.test.norm[,-1]
# Extrahiere X und Y Daten
X.train <- df.train.norm[,-1]
Y.train <- df.train.norm[,1]
X.test <- df.test.norm[,-1]
Y.test <- df.test.norm[,1]Build Model
Mit der Funktion fit_lstm kann ich nun das LSTM aufbauen. Hier habe ich zusätzlich zu der Epochengröße und der Batchsize noch viele weitere Hyperparameter. Man kann ein validaton_split, eine k.fold crossvalidierung, dropout layer, hidden layers, outputlayer oder den optimizer einstellen. Da es bei neuronalen Netzen keinen einzigen richtigen Weg gibt, muss man verschiedene Konstellationen ausprobieren, um die richtige Wahl der Hyperparameter zu finden. Bei dem Netz im unteren Codechunk habe ich mich für einen Mix von verschiedenen Hyperparametern entschieden und habe dadurch einen RMSE von 50 erreicht. Das ist eine durchschnittliche Fehlerquote von 50 Personen pro Tag, bei im Durchschnitt 3900 Personen pro Tag eine Abweichung von 1,28%. D.h., dass Netz macht seine Vorhersage mit 98,72 % richtig. Beim Trainieren ist mir aufgefallen, dass bei größerer Epochenanzahl auch bessere Ergebnisse zu erzielen sind. Ich denke, dass hier sogar noch bessere Ergebnisse möglich sind, wenn man mit mehr Rechenpower in der Cloud diese Netze laufen lässt.
# Build and fit LSTM model using k-fold cross validation rmse = 50 // 35 // 63
lstm <- fit_lstm(X = X.train, y = Y.train,
timesteps = ann.timesteps,
validation_split = 0.05,
epochs = ann.epochs,
batchsize = c(NA,ann.batchsize),
k.fold = 3, k.optimizer = "min",
hidden = data.frame(c(50,25),c("tanh")),
dropout = c(0.5,0.1),
output_activation = "linear",
stateful = F, return_sequences = T,
loss = "mean_squared_error",
optimizer = "Ftrl",
metrics = c('mean_absolute_error'))Speichere das Model.
Hier eine andere Variante, die zum Absturz von RStudio geführt hat.
Modellvergleich
Stündliche Vorhersage
In diesem letzten Schritt werden die RMSE Werte der linearen Regression und des LSTM Netzes miteinander verglichen.
# Vergleich RMSE
print(paste("Linear Regression RMSE =", round(hourly_linreg_test_rmse, digits = 2)))## [1] "Linear Regression RMSE = 152.69"
## [1] "LSTM RMSE = 4.25"
Es ist klar zu sehen, dass das künstliche neuronale Netz ein weitaus besseres Ergebnis liefert, als die lineare Regression. Dies entspricht den Anfangserwartungen, da bereits früh zu sehen war, dass viele Bezüge zwischen den Variablen nicht linearer Natur sind. Das neuronale Netz kann die Zusammenhänge anders bewerten und bezieht ebenfalls die Komponente der Zeit mit ein. Besonders deutlich wird dies im visuellen vergleich. Es ist klar zu sehen, dass die LSTM Vorhersage beinah deckungsgleich mit den Originaldaten verläuft, wohingegen die Vorhersagewerte der linearen Regression die Schwankungen zwar mitnehmen, jedoch viel zu wenig ausgeprägt.
ggplot() +
labs(title = "Vergleich Der Vorhersagen Für Stündliche Nutzerzahlen") +
xlab("Date") + ylab("User Count") +
geom_line(data = hourly_graphdata, aes(y = Value, x = Date, color = "black")) +
geom_line(data = hourly_graphdata, aes(y = Predictions, x = Date, color = "cadetblue")) +
geom_line(data = hourly_graphdata, aes(x = Date, y = LinReg, color = "orange")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))+
facet_wrap(month(hourly.df.test$start_datetime, label = T, abbr = F), scales = "free_x") +
scale_color_manual(name = "Legende", values = c("black"="black", "cadetblue" ="cadetblue", "orange" = "orange"),labels = c("Tatsächliche Nutzer", "LSTM Vorhersage", "Lineare Regression Vorhersage"))Tägliche Vorhersage
Im unteren Plot sieht man wie der Validation Score mit der Anzahl an Epochen abnimmt.
average_mae_history <- lstm$avg_qual
require(ggplot2)
ggplot(average_mae_history, aes(x = epoch, y = validation_qual)) +
labs(title = "LSTM: validation scores (MAE) per each epoch") +
xlab("Epoch") + ylab("Validation score") +
geom_line(color = "darkblue") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))ggplot(data = average_mae_history, aes(x = epoch, y = validation_qual)) +
labs(title = "LSTM: validation scores (MAE) per each epoch") +
xlab("Epoch") + ylab("Validation score") +
geom_smooth(color = "darkblue") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))Nun machen wir die Vorhersage und berechnen den RMSE. Außerdem denormalisieren und dedifferenzieren wir, damit wir die echten Werte haben zum Präsentieren.
# Nutzerdifferenzen vorhersagen
lstm.X.test <- as.LSTM.X(X.test, ann.timesteps)
predictiondiff <- lstm$model %>% predict(lstm.X.test, batch_size = ann.batchsize)
predictiondiff <- denormalize(predictiondiff[,1], nd$min[1], nd$max[1])
# Differenzierung rückgängig machen
testseries <- c(y_prior_first_test, y_test)
testseries <- testseries[-NROW(testseries)]
predictions <- invert_differencing(predictiondiff, testseries)
# RMSE
lstm_rmse <- rmse(testseries, predictions)
lstm_rmse## [1] 80.11411
Zuletzt werfen wir einen Blick auf den Plot. Dort sind die vorhergesagten Nutzerzahlen und die echten Nutzerzahlen abgebildet. Man sieht, dass das Netz sehr gut die Schwankungen abbildet und auch sehr genau die Werte Vorhersagt.
# Vorhersage für tägliche Nutzerzahlen für ein Neighborhood
graphdata <- data.frame(Date = df.test$datehour, Value = testseries, Predictions = predictions, LinReg = lin_reg_q2n_predictions[4:NROW(lin_reg_q2n_predictions)])
require(ggplot2)
ggplot() +
labs(title = "LSTM Forecasting for Bikeusers") +
xlab("Date") + ylab("Bikeuser") +
geom_line(data = graphdata, aes(y = Value, x = Date, color = "black")) +
geom_line(data = graphdata, aes(y = Predictions, x = Date, color = "cadetblue")) +
geom_line(data = graphdata, aes(x = Date, y = LinReg, color = "orange")) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black")) +
scale_color_manual(name = "Legende", values = c("black"="black", "cadetblue" ="cadetblue", "orange" = "orange"),labels = c("Tatsächliche Nutzer", "LSTM Vorhersage", "Lineare Regression Vorhersage"))Man sieht, dass die Graphen fast übereinander liegen und die Vorhersage mit den echten Nutzerzahlen fast identisch ist. Mit dieser Vorhersage könnte man nun die logistischen Operationen hinter einer Bikesharingflotte optimieren und so die Kundenzufriedenheit und Nutzbarkeit verbessern.
Ausblick
Die Lösung für einen Bezirk ist zufriedenstellend. Um eine Lösung für alle Bezirke oder sogar alle Stationen zu erstellen, gibt es mehrere Möglichkeiten. Die erste Möglichkeit ist, für jeden Bezirk ein eigenes LSTM Netz zu erstellen, das ähnlich dem hier Vorgestellten aufgebaut ist. Dies kann analog auf Stationsebene erfolgen, da so das Problem der zu Großen Datenmengen in einer Rechenoperation umgangen wird. Des Weiteren ist es möglich nur ein künstliches neuronales Netz zu erstellen, das zeitgleich eine Vorhersage für alle Bezirke trifft. Hierfür muss das Vorgehen geringfügig angepasst werden. Ist in dem hier vorgestellten Netz die Target-Variable ein Vektor, so muss diese in einem gemeinsamen Netz eine Matrix sein, bestehend den Nutzerzahlen aller Bezirke. So entstehen mehrere Output-Units in einem Netz. Diese Lösung erfordert allerdings mehr Rechenleistung, als ein LSTM für nur einen Bezirk. Theoretisch ist dies auch auf Stationsebene möglich. Dies würde eine enorme Rechenleistung erfordern. Um dennoch eine Umsetzung zu ermöglichen, kann es sinnvoll sein eine Cloudplattform einzubinden. So wird genügend Rechenleistung zur Verfügung gestellt, um die entstehenden Datenmengen zielorientiert zu verarbeiten und eine zeitnahe Lösung zu gewährleisten.
Abschließend lassen sich beide Forschungsfragen zufriedenstellend beantworten. Sowohl tägliche, als auch stündliche Nutzerzahlen können mit Hilfe eines künstlichen neuronalen Netzes zuverlässig vorausgesagt werden. Dies ist bisher nur für einen Bezirk erfolgt, lässt sich aber auf alle anderen Bezirke, sowie Stationen übertragen. So kann eine genaue Vorhersage der Nutzerzahlen zur besseren Ressourcenplanung genutzt und die Kundenzufriedenheit gesteigert werden.
Sources
https://www.kaggle.com/mathijs/weather-data-in-new-york-city-2016 https://www.citibikenyc.com/system-data https://www.kaggle.com/meinertsen/new-york-city-taxi-trip-hourly-weather-data https://www1.ncdc.noaa.gov/pub/data/cdo/documentation/LCD_documentation.pdf https://www.ncdc.noaa.gov/cdo-web/confirmation https://www.ncdc.noaa.gov/cdo-web/datasets/LCD/stations/WBAN:94728/detail https://www.eecs.tufts.edu/~dsculley/papers/ad-click-prediction.pdf https://keras.io/api/optimizers/ https://keras.io/api/optimizers/ftrl/ https://colah.github.io/posts/2015-08-Understanding-LSTMs/